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1. Introduction: 


Optimal trajectory problems usually involve the determination of a vehicle 
acceleration history that will accomplish a required change in vehicle 
state with a minimum expenditure of fuel. Functional optimization problems 
of this type can be reduced to boundary value problems in ordinary differen- 
tial equations by application of the well-known necessary conditions of 
optimality in the form of the classical Calculus of Variations, or one of its 
more recent counterparts, dynamic programming^ and the maximum principle 
However, except in rare instances, the principal parts of the solutions to 
these boundary value problems are not available in closed form. State-of-the- 
art guidance schemes ^ ^ ^ circumvent this difficulty by considering 
approximate formulations of the original problem that allow analytic construc- 
tion of major parts of the solution, so that only a simple iterative process is 
required. Since these approaches avoid some of the time-consuming numerical 
integration procedures that would be required to compute a general solution 
to the fundamental problem, the speed needed for real-time application has 
been the primary motivation for semi-explicit methods of this type. As a 
result of the approximations, the accuracy and flexibility of current flight 
schemes are limited, primarily in that they are nearly optimal only for short 
arcs of powered flight and for specialized mission (boundary-value) conditions 
in a restricted coordinate. system. This limitation can be relaxed somewhat 
in practice by use of special purpose adjustments, but only at the expense of 
additional preflight analysis and simulation. 


(1) DREYFUS, S.E. Dynamic Programming and the Calculus of Variations , 
Academic Press, Inc., New York and London, 1965. 

(2) PONTRIAGIN, L.S. et al. The Mathematical Theory of Optimal Processes , 
Interscience Publishers, John Wiley and Sons, Inc., New York and Lon- 
don, 1962. 

(3) SMITH, I.E. "A Three Dimensional Ascending Iterative Guidance Mode," 
NASA-MSFC Report MTP-AERO-63-49 , June, 1963. 

(4) MacPHERSON, D. "An Explicit Solution to the Powered Flight Dynamics 
of a Rocket Vehicle," Aerospace Corporation, Report TDR-169 (3126)TN-2, 
October, 1962. 

(5) CHERRY, G. W. "A General Explicit Optimizing Guidance Law for Rocket- 
Propelled Spaceflight," AIAA Paper 64-638, August, 1964. 
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General iterative procedures for solving two-point boundary value 
problems may be classified under two main headings: "direct" methods 

and "indirect" methods. Roughly speaking, direct methods search over 
the space of functions satisfying the boundary value requirements for a 
function satisfying the differential equations, whereas indirect methods 
search over the space of functions satisfying the differential equations 
for a function satisfying the boundary- value requirements. 

Prior to 1965, general (flexible) numerical procedures for computing precise 
optimal trajectories were far too unreliable in convergence and costly in 
computation requirements to be considered for real-time guidance. However, 
an indirect method for computing optimal trajectories, OPGUID, was de- 
veloped in 1965^’^ incorporating improved techniques to obtain a sub- 
stantial gain in speed, convergence, and flexibility. The principal 
improvements were the use of an efficient integration algorithm that was 
tailored to special features of the initial-value problem, and the use of 
closed form representations of the final-value transversal! ty conditions 
for general orbital injection missions of interest. 

Specifically, the OPGUID algorithm developed in 1965 required less than one- 
half second per iteration of the boundary-value problem (on an IBM 7094) as com- 
pared to a minute or more per iteration required by other existing algorithms 
that were applicable to launch vehicle trajectory problems at that time. A 
simple scaling rule for the amount of the Newtonian correction that was per- 
mitted per iteration resulted in an extremely large region of convergence that 
was surprisingly insensitive. to accurate initialization (e.g. a complete 180° 
reversal in the initial thrust direction could be corrected for in ten to twenty 
iterations for typical Saturn and Apollo launch trajectories). 

The indirect approach is particularly well suited for real-time use, where the 
guidance scheme must continually adjust to perturbations in initial conditions 


(6) BROWN, K.R. and JOHNSON, G.W. "Optimal Guidance for Orbital Transfer," 

IBM Report //65-221-0003H, Huntsville, Alabama, 30 August 1965. 

(7) BROWN , K.R. and JOHNSON, G.W. "Rapid Computation of Optimal Trajectories," 
IBM Journal of Research and Development , Volume 11, Number 4, July 1967, 
pp. 373-382. 
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which is readily accomplished by a single Newtonian iteration on the boundary- 

value problem. Feasibility of the use of OPGUID as a real-time guidance 

scheme for optimizing single-burn-arc (i.e. a sequence of thrusting stages 

separated by relatively short, fixed, staging intervals) orbital injection 

( 8 ) 

missions was demonstrated in 1966 v . 

However, many orbit transfer problems require the use of several burn arcs 
separated by relatively long optimal coast arcs. Several authors ^ 
have reported on the need to use the more complex "direct" methods, quasi- 
linearization or generalized Newton-Raphson, in order to obtain convergence 
for this problem for the restricted case of motion in a plane and fixed 

boundary conditions. However, a multi-burn- arc version of OPGUID, developed 

( 12 ) 

in 1967 , demonstrated that the attractive fundamental approach of OPGUID 

could successfully converge a general formulation of this problem, with variable 
boundary conditions. A sophisticated version of the multi-burn program (SWITCH) 
has been developed under Contract NAS 8-21315, that has successfully converged 
a variety of orbital transfer problems with an efficiency and reliability compar- 
able to that of the original OPGUID. Currently a maximum of 0.25 seconds compu- 
tation time is required per iteration on a CDC 6600 and normally only three to 
six iterations are needed for typical problems. 

As a result, the indirect method of SWITCH is not only feasible but considerably 
superior to existing implementations of -quasilinearization in convergence as 
well as efficiency. A principal feature of SWITCH is the use of classical two- 
body theory to render the computations for coast arcs explicit. Since high- 
thrust multi-burn orbit transfers usually involve coast arcs many times longer 
in duration than burn arcs, this results in a substantial saving in computation 


(8) BROWN, K.R. and JOHNSON, G.W. "Real-Time Optimal Guidance," IEEE Transac- 
tions on Automatic Control , Vol. AC-12, No. 5, October, 1967. 

(9) KENNETH, P. and McGILL, R. "Two-Point Boundary-Value Problem Techniques.," 
Advances in Control Systems , Vol. 3, C. T. Leendes (ed.), 1966. 

(10) McCUE, G.A. and BENDER, D.F. "Satellite Orbit Transfer Studies," NASA 
Report #N66-36050, 1966. 

(11) O'MAHONY, M.S., ESKRIDGE, C.D. and HANAFY, L.M. "The Optimal Solution of 
Trajectory Problems Consisting of Several Extremal Subarcs by the Gener- 
alized Newton-Raphson Algorithm" American Astronautical Society, Paper 
AAS 67-348, 1967. 

JOHNSON, G.W. and SHULL, N.W. "Optimal Guidance with Controllable Pro- 
pellant Mass Flow Rate," IBM Report CESD #009 December, 1967. 
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( 12 ) 


( 13 ) 

per iteration, A universal variable formulation of the two-body problem 
with closed-form expressions for the state transition matrix is used. This 
formulation was adapted in a novel way to avoid the cumbersome computation 
of the three dimensional tensor of second partial derivatives of final state 
with respect to initial state that is required when computing the partial 
derivative of final costate with respect to initial state. 

Since all the necessary partial derivatives are available at each iteration, 
as was the case for the original OPGUID, the SWITCH algorithm is appropriate 
for computing real-time corrections to in-flight perturbations. 

Unlike the OPGUID algorithm, the SWITCH algorithm does require reasonable 
initialization. That is, it is not possible with SWITCH as it was with 
OPGUID to misalign the thrust direction by 90 or 180 degrees and retain 
convergence. However, rough estimates of impulsive solutions have proved 
more than adequate as initialization in every trial case. Moreover, deforma- 
tion of the desired mission characteristics by very large amounts - for ex- 
ample a change in the eccentricity of the destination orbit from 0.0 to 0.5 - 
has been successful even when initialization values were unchanged. This is 
evidence that the SWITCH scheme would have a large safety margin in its ability 
to reconverge guidance solutions in response to worst-case real-time perturba- 
tions. However, it remains to be shown by vehicle flight simulation tests that 
the SWITCH algorithm possesses the speed and convergence properties that are nec- 
essary for real-time guidance. Such a demonstration would represent a signi- 
ficant advance since present real-time guidance schemes are not capable of 
revising an entire multiburn trajectory in response to in-flight perturbations, 
but can o r» 1 v modify one burn arc at a time. 


(13) GOODYEAR, W.H. "Completely General Closed-Form Solution for Coordinates 
and Partial Derivatives of the Two-Body Problem," The Astronomical Journ- 
al , Vol. 70, No. 3, April, 1965, pp. 189-192. 
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II 


2. Formulation of the Multi-burn Optimization Problem 

In this section, different formulations of the fuel optimi- 
zation problem for multiple burn trajectories are considered. 

It is shown that certain usual idealizing assumptions lead to 
an ill-posed optimization problem for which no solution exists, 
and several ways are discussed for avoiding such difficulties by 
more realistic problem statements. 

An Idealized Problem Statement 

The equations of space vehicle motion in a central gravitational 
field may be expressed as 


ur cm u 



where r is position, v is velocity, the direction of the unit 
vector u/|u| is control, c is the exhaust velocity of the rocket 
engine and m is the rate of change of vehicle mass due to pro- 
pellant expenditure, which also is part of control and can be 
chosen within the limits a < m < 0. 


The instantaneous rate of cost, L, is to be -m; i.e. we wish to 
minimize mass loss (or maximize final mass), so, letting state be 
x = (r,v,m) and costate be p = 3J^/9x = (q,s,w), the Hamiltonian is 


2.2 H = L 4 * x = — m 4- q^v + s*^ ( | j U | ) + wm 

4 I r 1 3 m |u| y 

which is minimized over thrust direction if and only if 


2.3 
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min I m T 

2.4 u H = (-1 + w + -1 — 1—) in + q v - — ^-r-(s r) 

T=T " |r | 3 

Letting S (the switching function) be (1 - w - c|s|/m), we have 
that 2.4 is minimized over m if m = a for S _< 0, and m = 0 for 
S > 0. So 


2.5 H° = SaU(-S) + q T v ±^(s T r) 

r 

where U is a unit step function, U(x) = 0 for x < 0, and 

U(x) = 1 for x j> 0. The costate equations are easily seen to be 


q = r (-3y | r | 5 r T s) + s(p|r| 3 ) 


2.6 s = -q 



m 


It immediately follows that the time derivative of the switching 
function is independent of m, or 


2.7 


S = - — | s 

m dt 


We have already incorporated two idealizing assumptions that are 

conventional: 

(i) Apart from thrust acceleration, motion is Keplerian; hence 
usually periodic. 

(ii) Thrust is truly proportional to mass rate; hence, mass loss 
is zero when thrust acceleration is zero. 

We now add two more assumptions: 

(iii) No terminal constraint is time-dependent. 

(iv) There is no limit on (or penalty for) the number of separate 
coasts or burns. 
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Taken together, assumptions (i) - (iv) preclude the possibility 
of globally optimal solution trajectories for most missions- We 
can see this by arguing from the necessary conditions of opti- 
mality- First, observe that assumption (iii) implies that the 
Hamiltonian H = 0. Second, observe that if we consider the vector 
p f = (q,s,w-l) , ; then altering the costate vector by multiplying 
p ? by any positive scalar k does not affect the resulting trajec- 
tory or any of the necessary conditions- Therefore, of the seven 
degrees of freedom available in choosing the costate vector p, 
only 5 degrees of freedom are usable- for selecting different opti- 
mal trajectories. Clearly, given an initial state x^, all these 
5 degrees of freedom plus the freedom of choice of terminal time, 
t^, are needed in order for the problem of reaching a prescribed 
position and velocity r,v to be even locally well-posed. We shall 
see, however, using our assumptions (i) - (iv) , that true optimality 
would impose still further requirements on the costate vector p so 
that fewer than 5 degrees of freedom are actually available in choice 
of costate. Suppose there is a fuel optimal trajectory for transfer 
from an original orbit to a given destination orbit. Let r,v,m be 
any state on this optimal trajectory at which the switching function 
S is negative. We suppose, by assumption (i) , that r,v define an 
elliptical (or circular) orbit. Let r 1 ,v* be any other position and 
velocity on that elliptical orbit. Starting from r f ,v’ ,m, an optimal 
maneuver for reaching the given destination orbit must be simply to 
coast around to r,v and then follow the original optimal maneuvers 
from that point on. For, if there were a cheaper maneuver, it could 
have been combined with the first part of the original optimal tra- 
jectory to r,v,m plus a coast to r'jv'jm, thus obtaining a better solu- 
tion to the original problem. Now this maneuver of coasting from r T , 
v 1 , m to r,v,m and following the old trajectory thereafter, being 
optimal, must itself satisfy the necessary conditions 2.1 - 2.6. It 
can be easily verified, however, that any such trajectory agreeing 
with the original trajectory during a burn and at least one switching 
point (point at which S = 0) , must agree exactly even in costate with 
the original trajectory (apart from a positive scalar multiplication 
of the vector (q,s,w-l)). This is a contradiction, because on. the 
new trajectory we must have S = 0 at r,v,m whereas we supposed that S 
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was negative on the original trajectory at that point. The 
only way out of the contradiction is to suppose either that the 
original trajectory was not, in fact, optimal after all or that 
S is never negative along an optimal trajectory. In the latter 
case, S must either be always positive (no burns at all), or by 
2.7, we must have S = 0 whenever S = 0. But this latter requirement, 
again by 2.7, clearly removes at least one degree of freedom from 
the available choices of costate. 

A more direct way of seeing the typical nonexistence of mass optimal 

trajectories under assumptions (i) - (iv) is available if we accept 

two lemmas. Lemma 1: given an optimal impulsive solution to an 

orbital transfer problem, there is no finite burn maneuver with as 

low a cost. Lemma 2: sufficiently small single- impulse maneuvers 

are approached arbitrarily closely in cost by the best single finite 

burn, single-coast maneuvers for achieving the same orbit change. 

(14) 

Both of these lemmas are easily verified from Robbins . We can 
easily see, since any impulse can be broken up into a large number of 
very small impulses, that lemma 2 implies that any one- impulse transfer 
between elliptical orbits can be matched arbitrarily closely in cost 
by a sufficiently large number of finite burns and coasts. This, in 
turn, together with lemma 1 show that the cost of an optimal impulsive 
maneuver must be a greatest lower bound of the set of all achievable 
finite burn costs, but not itself a member of that set. 

Realistic Problem Statements 

The main result of assumptions (i) - (iv) that leads to nonexistence 
of optical trajectories is that, at any point in a trajectory, one can 
insert a coast of arbitrary duration without increasing the cost of that 


(14) H.M. ROBBINS. "An Analytical Study of the Impulsive Approximation," 
AIAA Journal , Vol. 4, No. 8, pp. 1417-1423, August, 1966. 
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trajectory at all. If there were a fixed terminal time by which 
the maneuver had to be completed, or a penalty for restarting the 
rocket engine, or a nonzero rate of cost during coasts, then there 
would be a limit to how far the expedient of inserting coasts could 
be pushed. In practical orbital transfer missions, of course, there 
are either time limits or engine restart penalties or coast duration 
penalties. For example, existing rocket engines actually lose a 
small amount of mass per second, even during coasts, through venting 
of vaporized propellant. Also, astronauts and even on-board equip- 
ment such as flight computers and telemetry electronics, can function 
in space only for limited times. There is also a penalty for restart- 
ing and stopping the rocket engine. In fact, even assumption (i) is 
not precisely satisfied because perturbations due to earth oblateness, 
third body effects, solar wind, and other flight perturbations cause 
the motion to be not precisely periodic, although it is very nearly 
so. It would, however, be extremely difficult and also unsatisfactory 
to try to obtain a unique solution only on the strength of these slight 
departures from periodicity. 

There are three attractive alternatives for changing the problem 
definition so that a unique optimum will exist: a) to impose a termi- 

nal time constraint, b) to impose a limit on the number of engine 
restarts, or c) to incorporate a non-zero (and non-negligible) rate 
of cost during coasts. Any one of these alternatives would be a suffi- 
cient device and may be genuinely appropriate depending on the mission 
being analyzed. It is important to observe that the solutions given by 
alternative a) and the solutions given by alternative c) form the same 
class of trajectories. For consider a solution given by alternative a). 
The necessary conditions for optimality are as before except that the 
Hamiltonian H is not zero, although it is, of course, a constant along 
the trajectory. It can be argued that in this case the constant H must 
be negative, so let Lq (a positive number) be -H. If we now consider Lq 
to be a rate of cost that is independent of thrust, and remove the time 
constraint , we get the same solution but to a different problem: the 



problem of optimizing over free terminal time in the presence of 
an added constant rate of cost equal to L^. So, insofar as the 
mathematical nature of the solution is concerned, it makes no 
difference if we use alternative a) or alternative c) • Alternative 

b) , however, gives solutions that are not also solutions to problems 
involving alternatives a) or c) . In fact, we can see that solutions 
given by alternative b) must satisfy all of our original necessary 
conditions 2.1 - 2.6 plus the condition that He 0. For these necessary 
conditions can be deduced using only the assumption that each switch- 
ing time (that is, time at which the engine is started or stopped) 

is optimized relative to the rest of the trajectory and that the 
trajectory thrust direction history is optimized relative to the 
switching times. However, alternative b) can be made to approach 
alternative c) while holding the time penalty fixed by increasing 
thrust levels. Thus for high thrust missions, alternatives a) ,b) , and 

c) may turn out to be practically equivalent. 

Numerical Considerations 

It should be recognized that even when we resolve the basic problem 
of non-existence of global optima by changing the problem definition 
according to alternative a), b), or c) , some results of the original 
ill-conditioning are still important. We know now that since every 
trajectory that is fuel-optimal with respect to the choice of each of 
a fixed number of switching times satisfies the necessary conditions 
2.1 - 2.6, there must be an infinite progression of solutions to the 
boundary value problem all satisfying 2.1 - 2.6, and each one cheaper 
than the last. Let P]_,P 2 >** # a corres P on ^ing progression of normalized 
initial costate vectors (p^ = the result of dividing p^ by its own magni- 
tude). Since the space of normalized costate is closed and bounded, it 
follows that this space contains a limit point p^. Now consider the func- 
tion x(t) = g(t,p) defined as the result of propagating initial state 
forward in time to time t from initial state x^ at t^ using initial 
normalized costate = p. 
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It is clear that in the vicinity of p some functions of the 

oo 

trajectory (e.g. its duration) become arbitrarily sensitive to 
small variations in p. Such unbounded sensitivities would be 
a serious hazard to any iterative scheme for obtaining numerical 
solutions, whether to the problem definition given by alternative 

a) , b) or c). It is desirable to choose an approach under which 
the procedure searches over a set of independent variables such 
that either the dependent variables are never extremely sensitive 
to the independent variables and conversely or the nature of the 
search confines the independent variables to a "safe ' 1 region within 
which extreme sensitivities in either direction are avoided. It is 
evident, however, that no choice of the independent, and dependent, 
variables for the search can avoid extreme sensitivities unless the 
range of values of the variables is confined to a region which ex- 
cludes trajectories for which p is near p^. 

Alternative b) has the advantage that a simple and intuitively mean- 
ingful constraint of limiting the number of burns to two, or, in some 
cases, three, can be applied to each individual iterate in the search 
so that the range of values of the independent and dependent variables 
automatically excludes by a large margin the vicinity of 
trajectories corresponding to p near p^. There is also a distinct 
advantage for the logical structure of an implemented algorithm for 
numerical solution of problem b) in that the number of switching 
points is fixed, instead of varying from iterate to iterate as it 
would under alternatives a) and c) . Since both these advantages carry 
real weight from a practical point of view, we have adopted alternative 

b) as the preferred approach. However, if for any reason an implementa- 
tion of alternative a) or c) is desired, most of the current algorithm 
would be transferable without modification to the revised problem. 
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Thus the problem to be solved consists of finding a sequence of 
burns and coasts that achieves the desired orbit with minimum 
fuel subject to the constraint that the number of separate burns 
is limited to k, where k is usually two, sometimes three. In some 
cases, it is desired to permit an initial coast prior to the first 
burn as well as one coast after each burn. In certain other cases 
(such as circular to circular coplanar missions) the initial coast 
is not permitted since its length would be indeterminate under the 
necessary conditions . 


3 . Iterative Solution of the Boundary Value Problem 

In section 2, both dynamical necessary conditions and necessary 
boundary conditions were discussed for several possible multi-burn 
optimization problems. In this section, a method is described in de- 
tail for the numerical solution of one of these boundary value prob- 
lems. The particular problem chosen, called alternative b) in sec- 
tion 2, requires that trajectory which achieves the desired orbit (or 
orbital characteristics) with minimum fuel expenditure subject to a 
limit on the total number of separate burn arcs. 

The dynamical necessary conditions for this problem are equations 2.1, 

2.3, and 2.6. The boundary conditions are given at the left end by the 
initial position, velocity, and mass of the vehicle and at the right end 
by six conditions defining the desired characteristics of the destination 
orbit, which are discussed in section 5 for several important mission types. 
The intermediate point constraints are that the switching function S be 
zero at each interior switching time. In addition we require that the 
switching function be positive on coasting arcs and negative on burn 
arcs. This last requirement can be used to detect the desirability of 
adding further burn arcs. 

The sequence of burns and coasts of the desired trajectory is partially 
prescribed in advance as part of the problem statement. The number of 
separate burn arcs may be chosen, and the initial arc may be designated 
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as a burn or a coast. In a fundamental sense, a genuine orbit transfer 
problem which begins at a prescribed position and velocity, requires an 
initial coast arc in order to be fuel-optimal. But two important cases 
exist in which the initial coast must be suppressed. First, in the case 
of a circular to circular coplanar mission, the length of an initial coast 
would be completely arbitrary in the sense that all initial coast durations 
are equally compatible with fuel optimality. Second, there are cases in 
which mission geometry constraints preclude an optimal initial coast because 
it would result in a fall into the lower atmosphere or even beneath the 
earth’s surface. 


For purposes of specifying precisely the independent and dependent variables 
of the iterative scheme for solution of the boundary value problem, it is 
necessary to treat separately the case where an initial optimal coast is 
permitted and the case where the trajectory is required to begin with a 


burn arc. 


Let tg be the initial epoch for which the initial position. 


velocity, and mass are given: r(tg) = r Q> = V 0’ m ^ t o^ = m 0* Let n 

be the number of separate burn arcs (usually two or three). Then, for a 

trajectory beginning with a coast arc, the switching times are t = start 
th th 

of j burn and = end of j burn for j = l,...,n. Combining the 

constancy of H with the condition that S = 0, we have, at each switching point, 

T T 3 

by equation 2.5, the condition that q v - ps r/|r| = H. Recalling that this 

same condition is the principal transversality condition in the single-burn 

T T 

case, we will call q v - ps r/|r| 


the transversality variable Tv and so write 
Now this gives only about half 
of the actual necessary conditions at switching times, for it is easily seen 
from 2.1 and 2.6 that Tv(t) is identically constant along any coast arc, so 
if Tv(t 2 j) = H, it is automatically true that TvCt^^^) = H also. However, 
additional necessary conditions can be deduced from the fact that the switch- 
ing function S must be zero at the end of a coast as well as at the beginning; 
thus “ SC^j) = 0 • Since S=i-w-c|s|/m and since m and hence w 


the conditions as TvCt^) = H, i = l,...,2n 


are zero during coasts, this implies that l s ( t 2 j^q^ 


" |s(t 2 J | = 0. 


Thus we 

obtain the following set of 2n - 1 necessary conditions at intermediate switch- 
ing points: 


i. 
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(a) T v ( 2 j _i> = Tv(t 2 j)> J = 1 » •••> n 

(b) | s (*-2^+1^ I = I s ^ 2 j j = — 1 


The advantage of this set of conditions compared to the more straight- 
forward requirement that H(t^) = 0 and S(t^) = 0, i = l,...,2n - 1 is 
that it becomes possible to completely drop the costate variable w from 
all equations. For we can regard w_^ as defined by the requirement that 
S(t^) = 0 and then all necessary conditions are implied by (a) and 
(b). 


For missions in which the trajectory is required to begin with a burn arc, 
the switching conditions are altered in that the t^j and t 2 j+i now squire 
conditions appropriate to the beginning and ending of burn arcs respectively 
and in that there is no necessary condition applied at the beginning of the 
first burn. Thus we have 2n - 2 conditions: 


(a') Tv(t 2 ) = Tv(t 2j+1 ), j = 1, n - 1 

(b ’ ) | s (t 2 j ) I ~ l s < t 2 j_iH = °> J = !>•••> n - 1 


Observe that we have one more switching condition on an n-burn mission 
beginning with a coast than on an n-burn mission beginning with a burn. 

This corresponds to the fact that when an initial coast is permitted, its 
duration is an added variable for optimization. 

In summary, we can view the boundary value problems as follows: 

For a mission beginning with a coast arc: Find values for the six com- 
ponents of initial costate and the 2n switching times t ;L’-* , > t 2n suc ^ that 
the result of integrating 2.1 and 2.6 forward from t Q to t^^satisfies six 
right end mission conditions at t > the 2n - 1 intermediate switching con- 
ditions given by (a) and (b), and the condition | u Q | — 1. For a mission 
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beginning with a burn arc: Find values for the six components of initial 

costate and the 2n - 1 switching times i such that the results 

of integration 2.1 and 2.6 forward from t Q to satisfies six right 

end mission condition, the 2n - 2 intermediate switching point condition 
(a 1 ), (b f ) and * the condition I u 0 1 = 1. 

To solve numerically either of these problems, we apply the multi-dimen- 
sional Newton .method. Let the six components of initial costate and the 
2n switching times be combined into_ a vector y of dimension 2n + 6, and let 
the six constrained right end variables, the 2n - 1 intermediate switching 
variables of (a) and (b), and the variable |u 0 | be combined into a 
vector z of dimension 2n + 6. Then if z* is a vector of desired values for 
the components of z, the boundary value problem for missions with initial 

coasts becomes the problem of finding a vector zero of z* - z(y). A similar 

definition applies to the problem when an initial burn is required but with z 
and y having dimension 2n + 5 instead of 2n + 6. 

The Newtonian iterative solution scheme is familiar; starting from an initial 
guess y^, we obtain successive estimates 

3 - 1 y ±+i = y ± y [ z * - z(y i ) ] K 

where K is a scalar between zero and one. When K is one, the scheme is a 
pure Newton 1 s method; when K is less than one, the adjustment 
represents an attempt to obtain ~ z ± + K ( z * “ z i) ' This is often 

useful when z* - z ^ would be too large a change in z to correspond approxi- 
mately linearly to a change in y. 

In order to be able to implement 3.1, it is necessary to be able to compute 
the vector z(y) and the matrix 3 z/3y as functions of y. In the present case 
this means the ability to compute a trajectory and its first variations with 
respect to initial costate and the switching times t^. The trajectory is 
computed as a sequence of burn arcs and coast arcs, burn arcs being computed 
by numerical integration of 2.1, 2.3, and 2.6, and coast arcs being computed 
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explicitly by a method discussed in section 4. First variations of burn 
and coast arcs are computed concurrently with the arcs themselves. 

The method of integration used for the burn arcs is a fourth-order 
Runge-Kutta scheme suited for equations such as 2,1 and 2.6 (and the 
associated equations of first variation) which can be expressed as 
second order differential equations with first derivatives absent. In 
order to estimate and control truncation error by altering step sizes, 
a Richardson extrapolation method is used with a special combination of 
Runge-Kutta steps consisting of three steps of size h of integration of 

2.1 and 2.6 together with one overlapping step of size 3h of integration 
of 2.1, 2.6. and the associated equations of first variation. The result- 
ing scheme has proved very efficient and reliable and is described in 
full detail in reference 7. 

A special consideration that deserves mention in connection with the 
computation of 9 z/d y is the correct generation of those columns of the 
matrix that represent partial derivatives with respect to switching times. 
It is well known that if we have a vector system of differential equations 

3.2 a = f (a) 


the partial derivative of the solution a(t) with respect to any variable 3 
on which the initial value a q depends, satisfy the equations 


3.3 


d 3 q(t) 
dt 


a 3 


a f (a) 

a a 36 


3 a(t) 


and the initial condition, when t Q does not itself depend on 3, is 


3.4 


a a(t) 


3 a 


3 3 1 t=t Q 3 3 


The problem of computing 3 zfa y for purposes of the iteration 3.1 consists 
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inainly of determining the partial derivatives of such a vector a 
consisting of both state and costate, with respect to the components 
of y, and this is accomplished by solving 3.3 explicitly on coast 
arcs and numerically on burn arcs. But in the case of those compon- 
ents of y that are switching times t^, the initial condition is not 3.4 
as usual. Instead the condition is as follows: 


3.5 


9 a(t) | 

3 t t ' t=t ± 




aCt^) 


To appreciate this, it is sufficient to observe that since what is 
desired is the initial value of the (continuous) solution to a differen- 
tial equation (3.3), we really want 


3.6 



which, under the assumption that a = f^(a) before t^ and a = after t^, 

gives the following relation between a at time t shortly after t and a at 
time T shortly before t^ 


3.7 a(t) = a(T) + f ) (t ± - T) + f 2 (a(? 2 ))(t - t ± ) 


where C^eCTjt^) and . Differentiation of 3.7 with respect to ^ 

yields 

3.8 3 = f^cK^)) - f 2 (a(5 2 )) + 0(t ± - T) + 0(t - t ± ) 

1 

which in the limit as T t^- and t t^+ gives 3.5. 


4. Efficient Coast Arc Computations 

The computations required for burn arcs of the multi-burn SWITCH 
algorithm are essentially the same as those of the OPGUID algorithm, 
for single burn arc missions. The chief remaining computational 
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section of the algorithm is that which updates state, costate, and their 
partial derivatives along a coast arc. Considerable effort was ex- 
pended in planning these computations in order to take full advantage 
of the computational economies offered by the known efficient (semi-) 
explicit schemes for computing state and its partial derivatives along 
coast arcs. After extended analysis, it was found that costate and 
its partial derivatives with respect to initial costate and state 
could be computed easily using the corresponding computations for state 
and its partials. 

For each coast, we need values of r, v, u, and u (where u = -s) at 
the end of a coast given their values at the beginning and given the 
duration in time of the coast. We require in addition the partial 
derivatives of the final values with respect to the initial values since 
these partial derivatives form a necessary link in the chain of computa- 
tions leading ultimately to the partial derivatives of all the dependent 
variables of the boundary value search, described in section 3 , with 
respect to the independent variables. For purposes of the coast arc 
computations we can ignore m and w because their final values are simply 
equal to their initial values and neither m nor w appears in the differ- 
ential equations for r,v,u and u which apply during coasts: 



4.2 U = r(3p|r| 5 r T u) + u(-y|r| 3 ) 

It turns out to be convenient in what follows to define state x as (r,r) 
and (modified) costate y as (u,u), even though the costate that falls 
out directly from the calculus of variations or the maximum principle 
is p = (u,-u). The advantage of y over p is that it bears significantly 
simpler relations to state x. 

Completely general closed-form solutions are well-known for the propa- 
gation of state x according to 4.1 and for the corresponding computation 
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of the state transition matrix d x(t) /d x(t Q ) • A good formulation is 


Goodyear' s^ 13 ^ which bases all computations on a given initial value 

for state x and a time interval t - t . What we need beyond what 
o o 

Goodyear's or other usual formulations give are these additional six 
by six matrices: 9 y(t)$ y (t Q .) , 9 y(t)/5 x(t o ) , and 9 x(t)/^ yCt^) . Now, 

since x is independent of y along a coast, we can see immediately that 
the last of these matrices vanishes. 


4.3 


9 x(t)^ y(t Q ) = 0 


l 

) 


First we observe that if we write 4.2 as a differential equation for 


the whole (modified) costate vector y, we obtain 


4.4 


d 

__ Y (t) 
dt 



y(t) = 


9 x(t) 
9 x(t) 


y(t) 


Hence the differential equations for the matrices 9y(t)/9y(t o ) and 
9x(t)/9x(t ) are the same 


d 9 y(t) = 9 x(t) 9 y(t) 
dt 9y(t o ) 9 x(t) 9 y(t Q ) 


_d 9 x(t) = 9 £(t) 9 x(t) 

dt 9 x(t ) 9 x(t) 9 x(t ) 

o o 


Since both state and costate transition matrices must be initially 
equal to the 6x6 identity matrix 1^, we obtain 


4.6 


9 y(t) _ 9 x(t) 

3 Y(t„) " 9 x(t Q ) 


= 4>(t,t o ) 


A useful property of the matrix <|> ( t , t Q ) is the symplectic property 
which may be written: 


4.7 


<f>(t,t o ) T J<l>(t,t o ) = J 
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where 

4.8 J = 

This property may be derived by observing that 4.7 holds trivially 
at sinrt« <J>(t Q ,t o ) « 1 ^ and by writing out 

<f)(t,t o ) T j<}.(t,t o ) 

and noticing that it vanishes. The symplectic property is principally 
useful in obtaining <J)(t,t o ) ^ =<J>(t o ,t) for, by 4.7 and 4. 

4.9 4>(t,t o ) -1 = J T <j>(t,t o ) T J 

which is a simple rearrangement (with some changes of sign) of 
the elements of ( 1 5 1 ^) itself. 

Now since y is linear homogeneous in y, y(t) = 9y(t)^y(t Q ) y(t Q ), i.e. 

4.10 y ( t) = (j> (t , t Q )y (t Q ) 



The remaining matrix to be computed, 3 y(t)^ x(t Q ) is most directly 
expressable from 4.10; thus 


4.11 




3 (t o ) 



This equation would be very cumbersome to evaluate directly because 
of the 6 x 6 x 6 tensor of second partial derivatives 
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I 


9<Kt,t o ) ik ^ 3 3 x-j (t) = 3 3 x^ (t) 

5x j (t o } ^j^o 5 


But, by equation 4.12, 4.11 can be rewritten thus: 


4.13 


Rt . ) ,] = y 3 

l3x(to)J £ **^ 0 * 
ij 


3 x i (t) 

3 x. (t ) 
3 ° 


which can be interpreted as the linearized change d <f>(t,t o )^j in 
4> ( t , t Q ) due to a change in x q , dx Q = y(t Q ). This linearized change 
dcf) can be computed with very slight added labor when <)> itself is 
computed. 


This computation scheme has been implemented in a computer program which 
is much simpler than and about one half the size of the only known earlier 

program for explicit computation of state, costate, and their transition matrices 

i _ (15) 

along a coast arc 


5 . Right End Conditions for Various Orbital Missions 

Regardless of the kind of orbit transfer mission, the left end boundary 
conditions are given by the position, velocity, and mass of the vehicle 
at the initial epoch t Q , and the intermediate switching point conditions 
are as given in section 3. These conditions leave, in each case, six degrees 
of freedom remaining at the right end, which must be eliminated by necessary 
conditions that define the particular kind of orbit transfer desired. When 
fewer than six conditions define the desired character of the destination 
orbit, supplementary transversality conditions, requiring that the uncon- 
strained orbit parameters be chosen optimally, are added to make a total of 
six independent conditions. 

The most basic mission is defined by desired values for a complete set of 
five orbital constants; that is, five independent functions of position 


(15) PRESTON, E.L . "Computer Program for State Transition Matrices," 

Fall, 1967 Co-op Period, Mechanical Engineering Department, Purdue 
University . 
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and velocity that are constant in the absence of thrust. There cannot, 
of course, be any set of six independent orbital constants in this sense, 
for that would imply that position and. velocity are themselves constant 
during unpowered flight. Thus the five orbital constant mission is the 
most fully defined orbit transfer mission that does not involve the time 
origin of the destination orbit. Five orbital constants completely deter- 
mine the locus of the orbit, but not the absolute time at which any par- 
ticular position is reached. 


The most basic six constraint mission is rendezvous, which can be 
thought of as a mission in which five orbital constants plus the time 
origin of the orbit are specified. Thus not only the complete locus of 
the orbit but also the locus as a function of absolute time is prescribed. 

Also of interest are several missions in which fewer than five orbital 
constants are specified. While a full set of five orbital constants de- 
fines completely the plane of the orbit and the size, shape, and orienta- 
tion of the conic within that plane, it is possible with fewer than five 
orbital constants to prescribe only part of the geometry of the orbit and 
optimize over the remainder. 


From the calculus of variations or the maximum principle it is known that 
if the mission requirements are prescribed values for only k functions of 
final state 

5.1 § ^ (Xf ) = 0 i — 1 , • • • , k 

where the state vector x has n > k components, then optimality with 

respect to the n - k remaining degrees of freedom in final state requires that 

the final costate vector p^ lie the space spanned by the gradients of the 

constrained functions. This means that if n - k vectors a^x), i=k+l,...,n, 

can be found such that the a^ span the space orthogonal to the space spanned 

by the gradients 9 g i^ of the g^, then p f must be orthogonal to a^Cx), 

3 x 

i = k + l,...,n. In equations, we have 


5.2 


P f T a.(x) =0 i 


k + 1, 
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5.3 


a_. (x) = 0; i = j = k + 1,..., 


3 g ± (x) 
9 x 


In the case of the fundamental five orbital constant mission, k = 5 
and n = 6, so only one vector a^(x) is needed and it must be orthogonal to 
the gradient of every orbital constant. By definition, a function g^(x) 
is an orbital constant just in case 


5.4 


0 = dF [ s i (x) l 


> g-j/x) 

9 x 


where x is computed for unpowered flight. 


5.5 


r = v; v = - 


yr 


Therefore it is immediate that one vector a^ (x) that will satisfy 5.3 when 
the g^ are all orbital constants is 


5.6 




Thus the transversality condition 5.2 for this mission can be written as 


5.7 


T 

’ £ v f 


■>/ r £ - 0 


which will be recognized as the same as the switching condition which was 
applied at the end of an intermediate burn during the trajectory, now applied 
also to the end of the last burn. That is, 5.7 is equivalent to 


5.8 


Tv(t f ) = 0 


The five constrained functions g^(x) can be expressed in various equivalent 
ways, but probably the simplest formulation is this: Let h = rxv be the 

angular velocity vector and let e - -[-fr + *?} ^ the" eccentricity 

vector" whose magnitude is eccentricity and whose direction is the direction 
of pericenter. Then define g^(x) , for i = 1,...,5 as follows: 
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5.9 


gjL (x) = h 1 - h ± * 
g 2 ( x ) = h 2 - h 2 * 
g 3 (x) = h 3 - h 3 * 
g 4 (x) = e ± - 

g 5 ( x > = e 2 “ e 2* 

where "*" means "desired value of." 

It is easy to verify that these five functions g^(x) are independent and 
uniquely define the values of all orbital constants except in the isolated 
case where h^(x) = 0, which may be accommodated by reordering the coordinate 
axes. Accordingly, the six right end conditions for the basic five orbital 
constant missions have been implemented as 5.7 and 5.1 for g^(x) defined in 
5.9. 


Right end conditidtis for the rendezvous mission are easier to formulate but 
slightly harder to implement. We require that at the end of the last burn, 
at tp the position and velocity of the vehicle must match the position and 
velocity of a "target body." The position and velocity of the target, R(T) 
and V(t), must be given as functions of absolute time t, but the target body 
as such may be fictitious; as, for example, in the case of a synchronous 
orbit injection mission in which no actual body yet occupies the desired 
orbit. The six right end conditions are 


5.10 


g i(x f ,tf) = r^tj) - R 1 (t f ) = 0 


h + 3 (x f ,t f ) = v i (t f> ~ V ± (t f ) = °- 


>i = 1,2,3 


To implement 5.10 requires an algorithm for computing R(t f ) and V(t^) for 
variable t^. Since R(t), V(t) represent a free (unpowered) trajectory, this 
can be carried out using parts of the computations described in section 4, 
provided R (t Q ) and V(t Q ) are supplied for some initial epoch t Q . 
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6 


Useful four, three, and two-orbital constant missions are defined and 
transversality conditions answering to 5.3 derived for them on pp. 502 
and 503 of reference 8. 


Test Results 

The convergence properties of the multi-burn orbital transfer routine 
SWITCH were evaluated using basic Apollo type missions. The missions 
involved transfer from low earth orbits of about 100 nautical mile alti- 
tudes to orbits with altitudes ranging from 200 nautical miles to 19,300 
nautical miles. Both high and low thrust vehicles were used in the evalua- 
tion of the program. 

The first mission consisted of transfer from a near-circular orbit with an 
altitude of approximately 100 nautical miles to a circular orbit at a syn- 
chronous altitude of 19,300 nautical miles. The vehicle was assumed to have 
an initial mass of 1.26 x 10^ kgm, an exhaust velocity of 4.15 km/sec and a 
mass rate of 2.24 x 10^ kgm/sec. In order to use the SWITCH program one 
must estimate the independent variables of the boundary value search (i.e. 
the initial costate vector and the lengths of the burn and coast arcs) 
and fix the number of burn and coast arcs. The number of separate burn arcs 
was set at four (coast-burn-coast-burn) and the initial costate vector and 
the lengths of arcs were chosen to correspond to values used in evaluation 
of an earlier program. The SWITCH algorithm converged to the solution in one 
iteration. 

Since the initialization for this first mission was essentially equal to the 
converged solution (presented in Table 1), the final orbit was deformed into 
an ellipse to test the program more fully. The deformation was accomplished 
by choosing the eccentricity vector, e (the vector whose magnitude is equal to 
the eccentricity and which points in the direction of the pericenter) , to be 
parallel to the injection vector obtained in the first case and with the 
magnitude of the vector allowed to vary from .05 to .5. The initial orbit, 
vehicle parameters, and initial estimates of costate and arc times were fixed. 
The results are presented in Table 1 and as expected the number of iterations 
for convergence increased with increase in the eccentricity of the final orbit. 
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Of special interest is the case where the eccentricity was set at .5. 

The converged lengths of the second coast and burn arcs were changed 
by approximately 35 percent and a significant change was also required 
in the initial costate vector. The algorithm was able to overcome 
these difficulties and converged in five iterations. 

The right end boundary constraints were changed to allow solution of 
rendezvous missions and the same test cases were used to evaluate the 
altered algorithm. The solutions obtained for the five constraint 
missions were used to specify the radius vector, velocity vector and 
time necessary to describe the target vehicle. As indicated by Table 1 
the rendezvous missions generally required a larger number of iterations 
to converge. The sensitivity of the algorithm to changes in the time 
frame of the target vehicle was also investigated. The final orbit ec- 
centricity was fixed at .05 and the target vehicle was allowed to be both 
late and early in achieving the fixed position and velocity. The plus and minus 
200 second change in orbital epoch time required, in each case, about the 
same number of iterations for convergence as the nominal epoch case. The 
direction of the time change did not seem to affect the convergence pro- 
perties. 

The second test mission consisted of transfer from a 100 n.m. circular 
orbit to a 19,300 n.m. circular orbit which was rotated 44 degrees out of 
the plane. The vehicle characteristics were chosen to be the same as used 
in the first cases and the algorithm was able to converge to the solution 
in five iterations. The five constraint solution was again used to define 
the position of the target vehicle as a function of time and the rendezvous 
form of the algorithm converged in eight iterations . The results are pre- 
sented in Table 2, and, as shown, large changes were required in the initial 
costate vector. 

The third mission consisted of a burn-coast-burn transfer from a 96 n.m. 

circular orbit to a 196 n.m. circular orbit with a relatively low thrust 

vehicle. The vehicle was assumed to have a mass of 9248 kgm, a mass rate 

of 2.31 kgm/sec and an exhaust velocity of 2.16 km/sec. The initial costate 

vector was chosen such that u was in the direction of the vehicle velocity 

o 
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vector and that u was in the direction of the acceleration vector 
o 

of the coasting vehicle. The entire initial costate vector was scaled 
such that the magnitude of the u q vector was one. The length of the 
first burn arc was determined from the amount of velocity change required 
to deform the 96 n.m. orbit into a 96 n.m. x 196 n.m. orbit. The coast 
arc was chosen to be one half of the period of the elliptical orbit and 
the length of the second burn arc was chosen to gain the velocity necessary 
to circularize the orbit at apogee. The problem converged in two itera- 
tions and required five percent changes in the lengths of the burn arcs as 
well as a . 7 change in magnitude of the final u vector. 

The last mission tested required a transfer from 81 n.m. x 120 n.m. orbit 
to a coplanar 210 n.m. circular orbit. The vehicle characteristics con- 
sisted of an initial vehicle mass of 14,486 kgm, mass rate of 29.36 kgra/sec, 
and exhaust velocity of 3.07 km/sec. The mission was chosen to consist of 
four spearate arcs and the initialization of the costate vector and the 
switching times was accomplished by a rough pencil and paper estimation of 
the impulsive solution similar to that for the third mission. In this case 
the length of the initial coast arc was required and it was chosen to insure 
that the initial burn arc was centered around perigee. The converged solu- 
tion was within ten percent of the estimated arc lengths and required a .07 
change in the magnitude of the final u vector. This test mission required 
four iterations for convergence. 

All cases in which even a crude impulsive solution was used for initialization 
resulted in essentially immediate convergence. Only when an initialization 
appropriate to one mission was carried over unchanged to a significantly altered 
mission did the program require an appreciable number of iterations. 
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Destination Orbit Duration of Burn and Coast Arcs Magnitude of Number of 

Specifications Changes in u Iterations 


Type 

of 

Mission 

Eccentricity 

Minimum 
Altitude 
(n. m. ) 

1st Coast 
Arc 
(sec) 

1st Burn 
Arc 
(sec) 

2nd Coast 
Arc 
(sec) 

Znd Burn 
Arc 
(sec) 

| AuJ | Au f | 

to 

Converge 


Five 

0 . 

19, 300 

399.83 

255. 82 

18, 729. 54 

129. 11 

. 24E-3 .29E-3 

i 

Orbital 

. 05 

18, 300 

401. 16 

253. 47 

17, 572. 73 

134. 87 

. 80E-2 . 26E-2 

3 

Constant 

. 1 

17,200 

402. 47 

251. 14 

16, 541. 93 

140. 57 

. 16E-1 . 55E-2 

3 


. 5 

11, 700 

412.63 

232. 95 

11, 155. 42 

184. 04 

. 78E-1 . 27E- 1 

5 


R 

E 

N 

o 

0. 

19, 300 

399. 83 

255. 82 

18, 729. 51 

129. 11 

. 24E-3 .29E-3 

1 

N 

M 

I 

. 05 

18, 300 

401. 15 

253. 47 

17, 572. 72 

134.87 

. 81E-2 . 26E-2 

3 

D 

E 

N 

. 1 

17, 200 

402. 48 

251. 14 

16, 541. 81 

140. 57 

. 16E-1 . 55E-2 

5 

Z 

T T 

A 

L 

. 5 

11, 700 

412. 61 

232. 95 

11, 155. 54 

184. 04 

. 78E-1 . 27E- 1 

28 

V 

o 

Target 200 
sec LATE 

. 05 

18, 300 

379. 53 

253. 80 

17,827.69 

134. 62 

. 17E-0 .29E-1 

3 

u 

Target 200 

. 05 

18, 300 

423. 07 

253. 80 

17, 316. 21 

134,63 

. 16E-0 ,19E-1 

4 

s 

sec EARLY 











FIRST ITERATION SPECIFICATIONS 400. Z55. 82 18, 7Z7. 46 1Z9. 11 


Table 1 


Low Altitude to Coplanar Synchronous Missions 



Duration of Burn and Coast Arcs 


Magnitude of 
Changes in u 


Number of 


Type 

of 

Mission 

1st Coast 
Arc 
(sec) 

1st Burn 
Arc 
(sec) 

2nd Coast 
Arc 
(sec) 

2nd Burn 
Arc 
(sec) 

1 Au 0 l 

lAUpl 

to 

Converge 



| 


Five 

Orbital j 

Constant 

1211. 55 

255. 41 

18, 728. 65 

124. 99 

. 58E-1 

. 78E-0 

5 

RENDEZVOUS 

1211.56 

255.41 

18, 728. 56 

124. 99 

.59E-1 

. 78E-0 

8 


FIRST 

ITERATION 1113.69 255. 18 18, 729.19 118.80 

SPECIFICATIONS 


Table 2 


Low Altitude to Synchronous with 44° 
Plane Change 










7 Conclusions 


The development of the OPGUID algorithm in 1965 provided an 
efficient and reliable means for optimizing single-burn-arc 
transfer missions. The basic approach of the OPGUID scheme 
was to use an ordinary shooting method but with carefully 
chosen coordinate systems and numerical methods aimed at 
speed and range of convergence. The simulated use of OPGUID 
as a real-time guidance scheme demonstrated superior perform- 
ance to the existing IGM semi-explicit scheme, especially in those 
cases involving large perturbations from the planned normal mis- 
sion characteristics. The ability to recover smoothly in real- 
time from large perturbations is due to the well-behaved nature 
of the OPGUID formulation and its avoidance of artificial approxi- 
mations that limit the flexibility and accuracy of current guid- 
ance schemes. 

The OPGUID scheme handled single-burn missions well, including 
relatively short coast arcs of fixed duration -as in the case 
of a multiple stage boost flight. General multiple burn prob- 
lems, however, require long coast arcs (often about one half of 
an orbital period) and permit optimal choice of the lengths of 
both burns and coasts. A different algorithm was needed to deal 
efficiently with the very long coast arcs and to choose optimally 
the associated switching times. 

A sophisticated multi-burn optimization program, SWITCH, has been 
developed under Contract NAS 8-21315, and has converged a variety 
of orbit transfer problems with an efficiency and reliability ap- 
proaching that of OPGUID. Computing time per iteration is under 
one quarter second on a CDC 6600, and most missions require fewer 
than half a dozen iterations. However, it remains to be shown 
that the SWITCH algorithm possesses the speed and convergence pro- 
perties required for real-time guidance. Such a demonstration would 
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be a significant advance in the art since present guidance schemes 
cannot revise an entire multiburn trajectory in response to in- 
flight perturbations, but can only modify each single-burn arc 
separately. In addition to adapting the multiburn optimization 
program for real-time use, it is desirable to extend the boundary 
value formulations beyond that of the five constraint orbital con- 
stant and rendezvous missions that have already been solved. 
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Appendix: SWITCH Program 


I. Users Guide 

The MAIN routine reads the data necessary to specify the mission. 
The data is read as follows: 


Card No. 

Variable Names 

FORMAT 

1 

NCASE 

12 

2 

UK, AMASS, BMAX, CEXV, HMAX 

5F10.9 

3 

XOE 

6F10. 9 

4 

QOE 

6F10. 9 

5 

C 

7F1 0. 9 

6 

LEGMAX, IMAX 

212 

7 

ATP(1 , 1), ATP( 1,2) 

2F1 0. 9 

8 

0 

ATP(2, 1), ATP(2 , 2) 

2F1 0. 9 

• 

8+LEGMAX 

• • 

• • 

ATP(LEGMAX+1 , 1) ATP(LEGMAX+1 , 2) 

2F1 0. 9 


NCASE appearing on card 1 tells the program how many separate 

data cases are to follow. Cards Z - 8+LEGMAX completely define a 

case and must be repeated for each separate mission. Card Z contains 

a description of the vehicle and the gravitation attraction of the earth. 

The variables are defined as follows: UK — gravitational constant 

3 Z 

(km /sec ), AMASS — initial vehicle mass (kgm), BMAX — mass flow 
rate (kgm/sec), CEXV — exhaust velocity of the engines (km/sec) and 
HMAX — maximum allowable integration step size (sec). XOE is a six 
vector containing the initial position of the vehicle in its first three 
components and the initial velocity in its last three components. These 
vectors are specified in km and km/sec in an inertial geocentric 
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cartesian coordinate system. QOE is a six vector containing the initial 

T 

costate vector (u , u ) with the magnitude of u assumed to be unity, 
o o o 

The C vector for the basic five orbital constant mission contains the de- 
sired angular velocity vector h = r x v of the final orbit in its first three 
components and e e (first two components of the eccentricity vector of 

the final orbit e = - —■ — — — j as its last two components. LEGMAX 

L I r I A* 

is an integer specifying the number of burn and coast arcs and IMAX is the 
maximum allowable number of iterations for one mission. The ATP array 
defines the duration of the coast and burn arcs. ATP(1, 1) contains the 
starting time of the mission and ATP(1,2) is positive if the first arc is a 
coast and negative if the first arc is a burn. ATP(2, 1) defines the time 
at which the first arc is completed and the second arc started. The last 
arc is always assumed to be a burn and therefore ATP(LEGMAX,2) should 
always be negative. 

In order for the program to perform properly, the mass, mass rate 
and exhaust velocity should be such that the vehicle is able to carry out 
the required mission. If the mass is very large in comparison to the 
thrusting ability c£ the vehicle it will take very long duration burns to 
move the vehicle from its present orbit. If the mass rate is large and the 
usable mass small, the amount of possible burn time would be very small. 
A reasonableness check should be made before submitting any test case 
to the computer, for experience has shown that a very large percentage of 
unsatisfactory computer runs are due to gross errors in the input data 
which tender the missions essentially impossible. 

II. Program Organization 

The multi-burn optimal guidance package consists of nine subroutines 
SWITCH, COAST, AMULT, RKG031, RKSTEP, YDDRHS, BVEVAL, 
ADJUST, and SOLVE. The entire package was written in FORTRAN IV 
and was tested on a CDC 6600 data processing system. The algorithm. 



not including the two routines MAIN and OUT used for inputting and 
outputting data, required approximately 7300 core locations and required 
at most . 25 seconds for each iteration. 

The SWITCH subroutine is the executive routine for the program. 

It accepts the input data from MAIN and performs the necessary initiali- 
zations. SWITCH computes an entire trial trajectory for each iteration 
including computation of the necessary partial derivatives by means of 
one call on COAST or RKG031 for each coasting or burning arc of the 
trajectory. At the end of each arc control is returned to SWITCH which 
then calculates the additional partial derivatives associated with the end 
of the arc as well as the difference between the desired condition at the 
end of the arc and the actual condition. When all arcs have been com- 
pleted it calls the BVEVAL and ADJUST subroutines which calculate the 
allowable changes in the independent variables. At the end of the iteration 
SWITCH checks to see if a solution has been reached or if the total num- 
ber of iterations equals the maximum that is allowable. If either of these 
conditions is met SWITCH returns control to MAIN and the mission is 
assumed completed. 

The COAST routine propagates the trajectory along coast arcs. It 
accepts as inputs the present state and costate of the vehicle as well as 
the length of the proposed coasting arc. It calls AMULT to do the neces- 
sary matrix multiplications and outputs the state and costate at the end of 
the coasting arc as well as the partials of final state and costate with 
respect to initial state. 

The burn subarcs are propagated by subroutines RKG031, RKSTEP 
and YDDRHS. RKG031 acts as the executive routine for burn subarcs. 

It initializes the arrays needed for the other routines, controls the integra- 
tion step size and determines when the burn arc has been completed. When 
the final burn arc has been completed, the routine also determines the 

partials of state and costate with respect to the final time, T . 

Jb 



The RKSTEP routine performs a single Runge-Kutta fourth-order 
numerical integration step on the system of second-order matrix differ- 
ential equations which describe the motion of the thrusting vehicle. The 
integration scheme requires three separation evaluations of the right hand 
side of the differential equation and these are performed by the YD.DRHS 
routine. 

When the entire trial trajectory has been calculated, the BVEVAL 
routine is called. This subroutine computes the difference DC between 
the desired end conditions C and their actual value as determined by XF 
and QF. It also computes the matrix G of partial derivatives of C with 
respect to XF and QF and the matrix E of partial derivatives of C with 
respect to QO and the switching times. The DC vector is then added as 
an additional column of the E matrix and the ADJUST routine is called. 
This routine immediately calls the SOLVE routine which performs a 
Gauss -Jordan reduction of the E matrix and determines the adjustment 
necessary, on a linearized basis, to null out the DC vector. In order to 
insure that the assumption of linearity is not violated, a quantity CK is 
calculated which limits the adjustments to be performed on the indepen- 
dent variables of the boundary value search. The ADJUST routine then 
calculates the new QO and switching times and returns control to SWITCH, 


III. Variable Definitions 


Program Name 
XOE(I) 

QOE(I) 

C(I) 

ATP(I, 1) 
ATP(I, 2) 


Mathematical 

Symbol 

x 


o 



C 


Description 

Initial state vector 

Estimated initial costate vector 

Vector of desired end conditions 

Switching times 

Type of ARC (+ Coast, -Burn) 
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Mathematical 


Program Name 

Symbol 

Description 

AM 

M 

Mass at end of each leg 

AMASS 

M 

Initial mass 

BMAX 

M 

Mass flow rate 

CEXV 

V EX 

Exhaust velocity of the engines 

HMAX 

h 

Maximum allowable integration 
step size. 

LEGMAX 

- 

Number of burn and coast arcs 

IMAX 

- 

Maximum number of iterations 


• 

allowable 

QOl(I) 

% 

Initial costate vector, updated 


U 

each iteration 

UK 

V 

Gravitation constant 

QO(i) 

% 

Costate at start of arc 

QF(I) 

q F 

Costate at end of arc 

XO(I) 

X 

o 

State at start of arc 

XF(I) 

F 

State at end of arc 


^ r 

Partials of end conditions with 

E(I,J) 

o ^ 

respect to the independent 
variables 

3 q oV ' ‘ t F 

DC(I) 


Vector difference between 

(Also last column 


desired end conditions and 

of E matrix) 


attained end conditions. 

Z(I,J) 

3 x,q 

Partial of present state and 
costate with respect to the 
independent variables 

3 q 0 ,t r** t f 


3 x -p 3 q f 

Partial of state at end of 

PHI(I, J) 

3 x 0 " 3 q 0 

coast with respect to state 
at start of coast 
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Program Name 


Description 


DPHI(I, J) 

DXF(I) 

DRF 

DVF 

CK 


Matheira tical 
Symbol 



Partial of costate at end of 
coast with respect to state 
at start of coast 


Projected change in final state 


Magnitude of projected change 
in position 


Magnitude of projected change 
in velocity 


Fraction of the current Newton - 
Raphson adjustment accepted 

(0 < CK < 1) 
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CO 


IV - PROGRAM LISTING ‘AND SAMPLE OUTPUT 


PROGRAM MAI N ( I NPUT * OUTPUT » TAPE 5 = INPUT *TAPE6=OUTPUT ) 
COMMON /CCP IN J /UK* LEG t ATP (7 ,2 ) * AM , BMAX ♦ CEXV 
* DIMENSION C( 12) »Q0E(6) ,X0E(6) * 

1 FORMAT ( 7F10.9) 

READ (5*3) NCASE 
NC = 0 

5 NC=NC+1 

WRI TE ( 6 ,6 ) NC 

6 F OR MAT(1H1»13H CASE NUMBER .12*//) 

READ (5 *1 )UK, AMASS* BMAX, CEXV. HM AX 
READ( 5,1 ) XOE 
READ ( 5,1 )QOE 
READ ( 5 * 1 ) ( C ( I) *1 = 1,7) 

READ (5*3) LEGMA X » I MAX 
FORMA T ( 2 1 2 ) 

FORMA T ( 2F10 • 9.) 

L1=LEGMAX+1 

READ ( 5 *4) (ATP (I *1 ) ,ATP (I * 2 ) * I = 1 . * Ll ) 

CALL SWITCH ( AMASS .LEGMAX , I MAX * C , QO E * XOE » HM AX ) 
IF(NCASE.GT.NC) GO TO 5 
STOP 
END 


SUBROUTINE OUT ( X * Q » X F , QF , LEGMAX , NPATH ) 

COMMON /CCPINJ/UK*LEG,ATP ( 7,2 ) , AM , BMAX , CEXV 
COMMON/ WAD J/ A ( 8 ) ,DXF(6) 

COMMON /WCOAST /PSY, ALPHA, FT 

COMMON/ WSWI T/E ( 12 ,13 ) ,DC (12 ) , DRF , DVF , CK , EVT , KCOUNT , BURN T 
DIMENSION X(6) *H ( 3 ) *E1(3) .TIME (50,7) » DC1 ( 50 » 1 2 ) , DRF1 ( 50 ) ,DVF1 ( 50 ) 
DIMENSION KC 1(50), 0(6), XF(6),OF(6),CK 1(50) 

I F ( NPATH ) 1,2,3 

1 LFG1=LEG-1 

IF(LEG.LT.3) WR I T E ( 6 , 40 ) KCOUN T 
40 FORMAT ( 1H0.17HITERAT ION NUMBER ,I3*/1 

IF(LEG.EO.l) GO TO 2 
W R I T E ( 6 , 1 0 ) L EG 1 , X , 0 

10 FORMAT (20X.9HC0AST ARC , / 2 5X , 4H LEG= ♦ 

1 I 2 , 1 X 1 2HST ATE AT END , 1 X , 6 E 1 4 . 6 , / 30X , 1 4HC0ST ATE AT END , 1 X , 6 El 4 . 6 ) 
WRI TE ( 6,22 ) PSY, ALPHA, FT 

22 FORMA T ( 30X ,4HPS Y= ♦ El 4 * 6 » 5X , 6HALPHA= ♦ E 1 4 • 6 » 5X , 22HCALCUL ATED COAST T 

1 IMF* ,F14.6) 

2 H( 1 ) =X ( 2 >*X ( 6) -X ( 3 ) *X ( 5 ) 

H( 2 ) =X ( 3 ) *X ( 4) -X ( 1 ) *X (6 ) 

H( 3 >*X (1 )*X( 5) -X(2 )*X(4 ) 

HM= SORT ( H ( 1 ) **2 + H( 2 ) **2 + H ( 3 ) **2 ) 

RM=SORT (X ( 1 ) **2+X ( 2 ) **2+X ( 3 ) **2 ) 

El ( 1 ) =-(X ( 1 ) /RM+ ( H ( 2 ) *X ( 6)-H( 3 )*X ( 5) ) /UK ) 

El ( 2 ) =-(X (2 ) /RM+ ( H( 3 )*X ( 4 ) -H ( 1 )*X(6 ) ) /UK) 

El ( 3 ) =-(X (3) /RM+ (H( 1 )*X ( 5) -H ( 2 ) *X<4 ) ) /UK ) 

EM = SORT ( El ( 1 )**2+El (2)**2+Fl(3)**2) 
ENERGY=-UK/RM+.5*(X(4)**2+X(5)**2+X(6>**2) 

AAXIS=-UK/( 2.0*ENFRGY). 

PMJN= AAXI S* ( 1.0-EM) 

RM A X = AAX I S* ( 1.0+EM) 

D FR I OD= ( 6,2831853 )*SORT ( ABS ( AAXIS**3/UK ) ) 

WRITE (6,11 ) AAX I S , RM IN, R MAX, ENERGY ,PERI OD ,HM ,H , EM*E1»RM 

11 FORMAT (25X.15HSEM I MAJOR AX I S= , E 1 4 . 6 , 1 X5HRM I N= , E 1 4 . 6 , 1 X 5HRMAX= , E 14 

1 .6, 1X7HENERGY=»E14.6»/25X,7HPERI0D=,E 14.6,1 X5HHMAG=,E 14.6* 1X8HH VE 
2CTOR ,3E14.6,/25X5HEMAG=,E14.6, 1X8HE VECTOR , 3E14 . 6 ♦ 1 X5HRMAG= , E 14 . 6) 
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TF< NPATH) 5*6*6 

5 WRITE(6,12)LEG,XF,QE,AM 

12 FORMAT (20X»8 HBURN A RC , / 2 5X , 4HLEG= , I 2 , 1 X 1 2 HSTAT E AT END , 1 X 6E14 . 6 ♦ / 
130X.14HC0STATE AT END * 1 X6E1 4 . 6 ♦ / 30X , 19HMASS AT END OF LEG=,E14.6) 

TF( LEG.LT.L.EGMAX) RETURN 
KC=MOD (KCOUNT-1 ,50) 

KC= KC + 1 
DO 7 1=1, LFG 

7 TIME(KC,I)=ATP( 1+1,1 )-ATP( 1,1) 

RETURN 

6 L6=6+LEG 

DETE=1.0 
DO 4 1=1, L6 
4 DETE=DETE*E( I , I ) 

TIME(KC»7)=BURNT 

WRI TE( 6,13) TIME (<C,7) ,(TIME(KC,I) ,1=1, LEG) 

13 FORMAT ( 20X »16HT0TAL BURN T I ME= , E14 , 6 , 1 X9HARC T I MES , 1 X4E1 4 , 6 , /6 1 X 
12F14.6 ) 

WRITE(6,14)DC*DETE,(F(T,I),I=1,L6) 

14 FORMAT (25X,2HDC,1X6E14b 6»/28X,6E14«6,/25X,17HDETERM!NANT OF E=, 
1E14. 6, 1X13HDI AGONAL OF E , 4E14 . 6 , / 2 5X , 8 El 4. 6 ) 

LI =LEG+1 
L2=L1+1 

WRI TE ( 6, 1 5 ) DXF *DRF»DVF»CK»EVT , (A<I)*I=1,L2) 

15 FORMAT ( 25X,4HDXF , 6 El 4 . 6 , / 2 5 X , 4HDRF= , E 1 4 . 6 , 5H DVF= , E 1 4 , 6 , 4H CK=, 
1E14»6,5H EVT=*E14»6*/25X*10HCK=MIN OF ,7E13.5) 

L7=LEG+7 

WRITE (6,30) (E( I ,L7) ,1=1 ,L6) 

30 FORMAT(/25X66HCHANGE REOUISTED IN INITIAL COST ATE , SW I TCH I NG TIMES 
1 AND FINAL T I M E * /2 5 X , 6Ei 4 . 6 , /3 5X , 6 E 1 4 , 6 ) 

WRITE(6*16) KCOUNT ,OF , < AT P( 1,1), 1=2, LI) 

16 FORMAT ( IX , 24HEND OF ITERATION NUMBER , I 3 , / 1 5X ♦ 7HNEW 00 ,6E14«6»/ 

1 1 5X ♦ 16HNEW SWITCH T I MES , 1 X6 El 4 * 6 , / / ) 

DO 8 1=1,12 

8 DC1 (KC, I )=DC(l‘) 

CK 1 ( KC ) =CK 
DRF 1 ( KC ) =DRF 
DVF 1 (KC)=DVF 
KC1 (KC)=KCOUNT 
RFTURN 

3 IF ( KCOUNT, GT. 50 ) KCOUNT = 50 

WRITE (6,17) 

17 FORMAT(1H1,50X,1 4HSUMMAR Y TARLES » / / , 1 X9HITERAT ION , 1X1 5HTOTAL BURN 
1 TTME,21X,29HLENGTH OF BURN AND COAST ARCS , /2X6HNUMBER »// ) 

DO 9 1 = 1, KCOUNT 

9 WR1'TE( 6,18 ) KC1 ( I ) ,TIME( I ,7) , ( TIME( I , J) , J=1 ,LEG) 

18 FORMAT ( 3X* I3»5X,7E15,7) 

WRI TE( 6, 19) 

19 FORMAT ( 1H0,// 1X9HITERATION »39X , 36HERROR IN BOUNDARY COND I T I ONS-DC ( 
11-8) ,/2X»6HNUMBER,//) 

WRITE(6,20) (KC1( I ) ,(DC1 ( I ,J) ,J=1 ,8) , 1=1 , KCOUNT) 

20 F ORMA T(3XI3,5X,8E14,6) 

WRITE (6,21) 

21 FORMAT ( 1H0,// 1X9HI TFRATION ,26X ,8HDC ( 9-12 ) , 3 2 X , 3 HDR F , 1 2 X , 3 HDVF , 1 2 X 2 
1 HCK , / 2X6HNUMBER , / / ) 

WRI TE( 6,18 ) ( KC1 ( I ) , ( DC1 ( I , J) , J = 9,12 ) ,DRF1 ( I ) ,DVF1 ( I ) ,CK1 ( I ) , 1 = 1 , 

1 KCOUNT ) 

RETURN 

END 
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SUBROUTINE SWITCH! AMASS,LEGMAX , I MAX , C , OOE , XOE ,HMAX ) 

COMMON /CCP I NJ/UK, LEG* ATP ( 7.2 ) * AM , BMAX * CEXV 
COMMON/WSWI T/E < 1 2 » 13 ) »DC < 12 ) ♦ DRF ♦ DVF , CK » EVT ,KC0UNT , BURNT 
DIMENSION C(12) »X0 ( 6 ) »Q0 ( 6 ) »XF(6)*0E(6)*Z( 12*12) ♦X0E(6)»Q0 1(6) 
DIMENSION OOE (6) * DUMMY! 13 ) ,PHT(6,6>,DPHT!6,6>*DZ(12»12) 

KCOUNT = 0 
N0=0 

L7=7+LEGMAX 
L6=6+LEGMAX 
EVT = 1 . E-8 
KMAX= IMAX 
L1=LFGMAX+1 
DO 2 1=1,6 

2 O01U)=Q0E(T) 

WR I TE (6,1 00) UK, AMASS, BMAX, CEXV, HMAX 

100 FORMAT (1H0*1X»3HUK=»F10,2»1X*13HIN!TIAL MASS= , El 6 . 8 , 1 X , 1 OHMASS RAT 
1E=,E16.8,1X,17HEXHAUST VELOC I T,Y= , E 1 6 . 8 * 6H HMAX = F1 0, 3 , // ) 

WR I TE (6,121) IMAX ,'LEGMAX 
121 FORMAT (6H0IMAX=* I 3 , 5 X , 7HLEGMAX = , 13,/ ) 

WRITE (6,101 ) (ATP! I ,1 ) ,ATP< I ,2 ) , 1 = 1 ,L1 ) 

101 FORMAT! 1H0 , 20 X , AH T I M E , 2 0 X ,27HT Y PE OF ARC! + COAST,- BURN)/(20X, 
1E16.8,15X,E16.8) ) 

WRI TE ( 6,102 ) XOE »Q01*(C(I)»I = 1»7) 

102 FORMAT(17HOINITIAL STATE XO ,6F1 7. 7/16H ESTIMATED Q0»6F17.7/ 

1 1 6H0DES I RED FINAL C,7E16.8/> 

1 LEG= 1 

AM= AM ASS 
KCOUNT =K COUNT + 1 
BIIRNT = 0.0 
DO 3 1=1,6 
00(1) =001 ( I ) 

3 XO ( Ii =XOF ( I ) 

DO 4 1=1,12 
E ( I ,13 >=0.0 
DO 4 J= 1 ,12 
E ( I , J) =0.0 

4 Z ( I ,J>=0.0 
DO 5 1=1,6 

5 Z( 1+6,1 ) = 1 . 0 

IF ( ATP! 1 ,2 )) 6,7,7 

7 CALL COAST ( X 0 » 00 » XF ,0F , PH I ,DPHI » UK , LEG , ATP , NO ) 

L FG6 = 6 + LFG 
LEG5=LEG+5 
DO 30 1=1,6 
DO 30 J=1 ,LEG5 
DZ ( I , J)=0.0 
DZ ( I +6 , J ) =0 • 0 
DO 30 K = 1 * 6 

DZ ( I , J ) =DZ ( I , J ) +PH I ( I ,K) *Z( K, J ) 

30 DZ ( 1+6, J) =DZ( 1+6 , JJ+PHI ! I ,K )*Z ( K+6 , J)+DPHI { I ,K ) *2 (K,J ) 

DO 31 1=1,12 

DO 31 J=1 ,LEG5 

31 Z( I ,J)=DZ(I,J) 

UM= SORT ( OF ( 1 )**2+QF( 2 )**2+0F ( 3 ) **2 ) 

CBMU= (CEXV*BMAX) / ( AM*UM ) 

DO 8 1=1,3 

8 Z< I+3,6+LEG)=-CBMU*QF( I ) 

IF! LEG.LE.l ) GO TO 12 
E(6+LEG,7+LEGMAX) =UMP-UM 
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DO 9 J=1 »LEG6 
E ( LEG6, J) =-E (LEG6* J ) 

DO 9 K = 1 ♦ 3 

9 E(LEG6*J)=E(LEG6*J)+(QF(K) /UM ) *Z ( K+6 » J ) 

DO 11 1=1*6 
0° ( I )=QF( I ) 

11 XOII)sXF(I) 

LEG=LEG+1 

6 CALL RKG031 (XO*QO»XF»QF»Z»EVT »HMAX»LEGMAX»NO) 

BURNT=BURNT+ATP(LEG+1 *1 ) — ATP (LEG»1 ) 

CALL OUT (X0*O0*XF»QF ♦ LEG MAX * -1 ) 

LEG6=6+LEG 

IF( LEG. GE. LEGMAX ) GO TO 18 
13 UM= SORT ( OF ( 1 ) **2+QF ( 2 ) **2+QF ( 3 ) **2 ) 

(jMP = UM 

CRMU = ( CEXV*BMA X) / ( AM*UM ) 

DO 14 1=1,3 

14 7. ( f+3 , 6+LEG) =CBMU*QF(I) 

DO 15 J=1 »LEG6 

DO 15 K = 1 » 3 

15 E ( 7 + LEG, J ) =E ( 7 + LEG , J ) + ( QF ( K ) /UM ) * Z ( K+6 » J ) 

» 12 R2=XF(1 ) * X F ( 1 ) + X F ( 2 ) * X F ( 2) +X-F ( 3 ) * X F ( 3 ) 

RS=XF ( 1 )*QF ( 1 ) +XF ( 2 ) *QF ( 2 )+XF ( 3 ) *QF (3 ) 

C3=-UK/(R2*SORT(R2) ) 

C4=-3.0*C3*RS/R2 

E(6 + LFG*7+LEGMAX)=XF (4)*QF (4)+XF ( 5 ) *QF ( 5 ) +XF < 6 ) *QF ( 6 ) -C 3*RS 
DO 16 1=1*3 

ryjMMYf I )=-C3*QF( I )-C4*X (r ( I ) 

DUMMY ( 1 + 3 )=QF< 1 + 3 ) 

DUMMY ( 1+6 )=-C3*XF( I ) 

16 DUMMY ( 1+9 ) =XF ( 1+3) 

DO 17 J=1 *LEG6 

DO 17 K= 1 , 1 2 

17 E(6+LEG*J)=-DUMMY(K) *Z ( K * J ) +E ( 6+LEG , J ) 

DO 10 1=1,6 
00(1) =QF( I) 

19 XO ( I ) =XF ( I ) 

LFG=LFG+1 

IF(ATP(LEG,2 ) .GT.O) GO TO 7 

F(7,?)=0.0 

GO TO 6 

18 CALL RVEVALf XF,OF,Z *C,F*DC) 

DO 20 1=1*6 
16=1+6 

E ( I ,^ + LEGMAX )=DC( I) 

DC ( I 6 ) = E ( 1 6 » L EGMA X + 7 ) 

20 E( 6+LEG MAX* I ) =Q01 ( I ) 

CALL ADJUST (001 , E , Z * DRF , D VF * C K , LEGMAX * ATP ) 

EVT = AMAX1 ( 1 . E-14, 1. E-8*DVF**2 ) 

EVT =AM INI (EVT* 1 .E-8 ) 

CALL OUT(XF, OF, XO, 001 , LEGMAX, 0) 

TF( KCOUNT.GE.KMAX) GO TO 79 
IF ( DRF.GE..1 .OR.DVF.GE. .005 ) GO TO 1 
79 CALL OUT ( XO , QOE » XF , OF , LEGMAX * 1 ) 

RETURN 

END 


SUBROUTINE COA ST ( XO *Q0 , XF , OF * P H I * DPH I , UK ♦ LEG ♦ A TP ,N0 ) 
COMMON/WCOAST/PSY, ALPHA ,FT 
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COMMON /CCOA ST /ANN ( 2 * 2 ) *BNN (3*3) ,XX0(6«6)»R0<3),V0<3>,R<3)*V(3) 

1 ,DRO (3 ) *DVO(3 ) ,DR (3 ) *DV( 3) * DAN ( 2,2) ,DBNN< 3,3 ) *DXXO< 6,6 ) 

COMPLEX CALPH,CAPSY 

DIMENSION X0(6) ,Q0(6),XF(6) ,OF(6) .PHI (6*6) ,DPHI ( 6 , 6 ) » A TP ( 7 * 2 ) 
DIMENSION HO ( 3 ) 

LFG5=5+LEG 

C TIME OF COAST 

T=A TP ( LEG+1 » 1 )— ATP ( LEG »1 ) 

C 00 (Q AT START OF COAST) 

DO 40 1 = 1 ,3 
13=1+3 
Rom=xom 
vo< r ) = xo< r+3 ) 

DRO ( I )=Q0< I ) 

40 DVD ( I ) =00 (13) 

C RO , VO STATE AT START OF COAST 

JUMP = 0 

RMO =SORT ( RO ( 1)*R0(1 )+R0(2)* RO ( 2 ) +R0 ( 3 ) *R0 ( 3 ) ) 

DRMO = ( RO ( 1 ) *DRO ( 1 )+R0( 2 ) *DRO( 2 )+R0 ( 3 >*DRO< 3 ) ) /RMO 
SIGO = RO ( 1 ) *V0 ( 1 )+R0( 2 ) *V0 (2 )+R0 ( 3 )*V0< 3 ) 

DSIG0= (VO ( 1 ) *DRO < 1 )+V0 (2 )*DR0 ( 2 )+V0 ( 3 ) *DR0 ( 3 )+R0 ( 1 )*DVO< 1 ) 

1+RO ( 2 ) *DVO ( 2 )+R0 ( 3 )*DVO( 3 ) ) 

ALPHA=VO ( 1 )*V0( 1 )+V0<2 ) * VO < 2 ) +V0 ( 3 ) *V0 ( 3 ) -2 .*UK / RMO 
HO(l)=R0(2)*V0(3)-R0(3)*V0(2) 

HO ( 2 ) = R0 ( 3 ) *V0 ( 1 )-R0 ( 1 )*V0 ( 3 ) 

HO ( 3 ) =R0( 1 )*V0( 2 )-R0 ( 2 )*V0( 1 ) 

P0= (HO ( 1 ) *H0( 1 ) +H0 ( 2 ) *H0 ( 2 ) +H0 ( 3 )*H0(3 ) ) /UK 

80 RSY=T/P0 

I F ( ALPHA)81»82*82 

81 CALPH=CMPLX(SQRT(-ALPHA) *0. ) 

GO TO 83 

82 CALPH=CMPLX (0. *SORT( ALPHA) ) 

83 CAPSY=PSY*CALPH 
SO=REAL(CCOS (CAPSY) ) 

S1=REAL(CSIN (CAPSY) /CAL PH) 

S2= (SO-1 .0) /ALPHA 

S3 = (S1-PSY)/ALPHA 
FT=RM0*S1+SI GO*S 2+UK*S3 
RM=RM0*$0+SIG0*S1 +UK*S 2 
I F ( JUMP.FO.l ) GO TO 84 
PSY=PSY+(T-FT) /RM 
TF(ABStT-FT) .GE..0001) GO TO 83 
JUMP= 1 
GO TO 83 

84 FM1=-UK*S2/RM0 
F=1 • 0+FM1 

FD=-UK*S1/(RM*RM0 ) 

G=FT-UK*S3 
GDM1=-UK*S2/RM 
GD= 1 . 0+GDM1 
UKR3=UK/(RM*RM*RM) 

UKR03=UK/ (RM0*RM0*RMO ) 

DALPH=2.0*(VO( 1 )*DVO( 1 ) + VO < 2 )'*D VO < 2 ) +V0 ( 3 ) *DV0 ( 3 ) +UKR03* ( 

1 RO ( 1 ) *DRO ( 1 ) +R0 ( 2 )*DR0(2 )+R0( 3 ) *DRO (3) ) ) 

DAPA=DALPH/ALPHA 

DAPA2=DAPA/ALPHA 

DPSY=-(DRM0*S1+DSIG0*S2+RM0*(PSY*S0-S1 ) *DAPA* . 5 +S I GO* ( PS Y*S 1 *. 5- 
1S2 ) +DAPA+UK* ( PSY-1 . 5*S1+PSY*S0* . 5 ) *DAPA2 ) /RM 
DSO= ( ALPHA*DPSY+.5*PSY*DALPH)*S1 
DS1 = SO+DPS Y+ ( PSY+S0-S1 )*DAPA*. 5 
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DS 2 = S1 *DPSY+ ( • 5*PSY*S1— S2 ) *DAPA 
DS3=S2*DDSY+ ( PSY-1 .5*S1+.5*PSY*S0 )*DAPA2 
S4 = ( S2— PSY*PSY* • 5 ) /ALPHA 

DS4=S3*DPSY+( PS Y*PSY*. 5-2.0 *S2+.5*PSY*S1 )*DAPA2 
S5=( S3-PSY*PSY*PSY/6.0) /ALPHA 

DS5 = S4*DPSY+(PSY*PSY*PSY/6.0+< 2 .0*PSY-2. 5*S 1+. 5 *PSY*S0 ) /ALPHA) 
1 *DAP A2 

U=S2*FT+UK*(PSY*S4-3.0*S5 ) 

DU = DS2*FT+'JK*(DPSY*S4+PSY*DS4-3.0*DS5 ) 

DRM=SO*DRM0+DS0*RM0+S 1 *DSTG0+DS1*S IG0+UK*DS2 
DF= (-UK*DS2-FM1*DRM0) /RMO 
DG=-UK*DS3 

DGD=( -UK*DS2-GDM1*DRM)/RM 

roi=rmo*rm 

DR01=RM»DRMO+DRM*RM0 
DFD= (— UK*DS1— FD*DR01 )/R01 
DO 4 1=1*3 

DR ( I) =R0 ( I )*DF+VO ( I )*DG+DRO t I ) *F+DVO< I )*G 
OF ( I ) =DR ( I) 

R ( I ) = RO ( I )*F+VO( I ) *G 

DV( T ) = RO ( I ) *DFD+VO ( I ) *DGD + DRO ( I ) *FD+DVO ( I ) *GD 
OF ( T+3 ) =DV( I ) 

4 V( I ) = RO ( I ) *FD+ VO ( I ) *GD 

C R * V STATE AT END OF COAST 

IF(NO.EO.l) GO TO 90 

DUKR3=-3«0*IJKR3*DRM/RM 

DUR03=-3,0*UKR03*DRM0/RM0 

S1R0=S1/RM0 

DS1 R0=(DS1-S1R0*DRM0 ) /RMO 
S1R=S1/RM 

DS 1 R= ( DS1-S1R*DRM) /RM 
R02=l»0/( RMO* RMO ) 

R 2= 1 • 0 / ( RM*RM ) 

UUK03=-U*UKR03 
D(JUK3 = -DU*UKR03-U*DUR03 
ANN ( 1 , 1 )=-FD*Sl RO— FM 1 *R0 2 
ANN (1*2) =-FD*S 2 
ANN ( 2 * 1 J = FM1*S1 RO+UUK03 
ANN (2,2)=FM1*S2 
DUMMl=ANN ( 1*1 ) 

DAN (1*1) =-FD*DSlR0-DFD*SlR0+UKR03*DS2+DUR03*S2 
DUMMY=DAN( 1 , 1 ) 

DAN ( 1 ,2)=-DFD*S2-FD*DS2 
DAN ( 2 *1 ) =FM1*DS1R0+DF*S1R0+DUUK3 
DAN ( 2*2 ) =FM1*DS2+DF*S2 
CALL AMULT 
DO 5 1=1,3 
DO 5 J = 1 * 3 
DXXO ( I , J 1 =DBNN ( I *J) 

5 XXO ( I * J ) *RNN ( I ♦ J ) 

DO 6 1=1*3 

DXX 0 ( I ♦ I )=DXXO( I ,1 )+DF 

6 XXO ( I , I ) =XXO ( I * I )+F 
ANN(1,1)=ANN(1*2) 

ANN(2,1)=ANN(2»2) 

ANN (1,2) =-GDMl*S2 
ANN ( 2 , 2 ) =G*S2-U 
DAN (1*1 ) =DAN (1*2) 

DAN(2»1)=DAN(2*2) 

DAN ( 1 ,2)=-GDMl*DS2-DGD*S2 
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DAN ( 2 ,2 ) =-DU+DG*S2+G*DS2 
CALL amult 
DO 7 1=1*3 
DO 7 J= 1 » 3 
J3 = J+3 

DXXO ( I * J3 ) = DBNN ( I ,J ) 

7 XXO( I * J 3 ) = BN N t I *J> 

DO 8 1=1*3 
13 = 1+3 

DXXCM I * I 3 ) =DXXO ( I *I3)+DG 

8 XXO( I*I3)=XX0(I ,I3)+G 
ANN (2,1 )=-ANN( 1,1 ) 

ANN ( 2 , 2 ) =— ANN ( 1 , 2 ) 

ANN ( 1 *1. )=-FD*SlR-GDMl*R2 
ANN ( 1 ,2 ) =U*UKR3-GDM1*S1R 
DAN ( 2*1 )=-DAN( 1,1) 

DAN (2,2) =-DAN (1,2) 

DAN (1,1 ) =— DFD*S1 R-FD*DS 1R+UKR3*DS2+DUKR3*$2 

DAN ( 1 ,2 ) =-GDMl*DSlR-DGD*SlR+DU*UKR3+U*DUKR3 

CALL AMULT 

DO 9 1=1,3 

DO 9 J = 1 , 3 

13=1+3 

J3= J+3 

D X X 0 ( I 3 * J3 ) =DBNN ( I, J) 

9 XXO( I3-,J3)=BNN(T *J) 

DO 10 1=4,6 

DXXO ( I ♦ I) =DX XO (1,1) +DGD 

10 XXO( I , I ) = X XO ( I , I ) +GD 

ANN (1,2) = ANN (1,1) 

ANN (2,2 ) = ANN (2,1) 

ANN ( 2 , 1 ) =-DUMM 1 

ANN( 1,1 ) =-FD* ( S0/R01+R2+R02 ) -UUK03*UKR3 
DAN ( 1,2) =DAN (1*1) 

DAN (2,2) =DAN (2,1) 

DA N ( 2 ♦ 1 ) =-DUMMV 

DAN ( 1,1) =-UUK03*DUKR3-DUUK3*UKR3-DFD*( S0/R01+R2+P02)-FD*( ( DSO-SO* 
1DR01/R01) /ROl-2. 0* ( R 2 *DRM /RM+R 0 2 *DpM0 /RMO ) ) 

CALL AMULT 
DO 11 1=1,3 
DO 11 J = 1 , 3 
13=1+3 

DX X 0 ( 1 3 , J ) =DBNN ( I ,J) 

11 XX0( I3,J)=BNN( I, J) 

DO 12 1=1,3 
13=1+3 

DXXO( 13*1 )=DXXO( 13,1 >+DFD 

12 XXO( T3,I )=XXO( 13,1 )+FD 

C XXO PARTIAL OF XF WITH RESPECT TO XO 

C DXXO PARTIAL OF OF WITH RESPECT TO XO 

DO 50 1=1,6 

DO 50 J=1 ,6 

PHI ( I ,J) =XXO ( I , J) 

50 DPHI ( I ,J)=DXXO( I , J) 

90 DO 52 1=1,3 

XF ( I ) =R ( I ) 

52 XF ( I +3 ) =V ( I ) 

RETURN 

END 


A-13 



SUBROUTINE AMULT 

COMMON /CCOA ST / ANN ( 2 * 2 ) > BNN ( 3 , 3 ) , XXO ( 6 » 6 > » RO ( 3 ) » VO ( 3 ) » R t 3 > » V ( 3 > 

1 ♦ DR 0 ( 3 ) *DV0(3 ) »DR(3) *OV(3) »DAN(2*2) »DBNN(3,3) »DXXO<6*6) 

DIMENSION DA (3,2) 

DIMENSION A ( 3 * 2 ) 

DO 1 1=1,3 
DO 1 J = 1 * 2 

DA ( I * J ) =DR ( I ) *ANN ( 1 » J ) +DV ( I)*ANN< 2, J)+R ( I ) *DAN ( 1 ,J)+V( I )*DAN<2» J) 

1 A (I , J)=ANN( 1, J)*R ( I )+ANN (2 , J)*V( I ) 

DO 2 1=1,3 

DO 2 J = 1 ,3 

DBNN < T , J ) =A ( I ,1 )*DRO ( J)+A( I , 2 ) *DVO ( J ) +DA ( I , 1 ) *RO ( J ) +DA ( I ,2 ) *VO ( J ) 

2 BNN ( I , J)=A ( I , 1 ) *R0 ( J !+A( I » 2 ) *V0 ( J ) 

RETURN 

END 


SUBROUTINE RKG031 ( XO ,00 , X F , OF » Z , EVT , HM AX , L EGMAX , NO ) 

COMMON /CCPINJ/UK, LEG, ATP (7,2 ) , AM , 8MAX , CEXV 

DIMENSION X 0 ( 6 ) *Q0(6)*XF(6) »QF ( 6 ) ,2(12,12) » YN ( 6 , 1 3 ) ,YDN( 6, 13) , 
1 Y3H ( 6 * 1 3 ) ,YD3H(6,13) * YM ( 6 * 1 3 ) , Y DM ( 6 , 1 3 ) ,EY(6) ,EYD(6) 

LEG6=LEG+6 
DO 1 1=1,3 

YN (1*1 )=X0( I ) 

YN ( 1+3 ,1 ) =00 ( I ) 

YDN ( I ,1)=X0( 1 + 3) 

YDN ( 1+3,1 ) =Q0 ( 1+3) 

DO 1 J=2,LEG6 
YN( I , J) =Z ( I , J— 1 ) 

YN ( I +3 » J ) =Z ( I +6 , J-l ) 

YDN ( I , J) =Z( 1+3 , J-l ) 

1 YDN ( I+3,J)=Z ( I+R* J-l ) 

TF=ATP(LEG+1 ,1 ) 

T0= ATP ( LEG, 1 ) 

H=SIGN(HMAX,TF-TO) 

TN = TO 
GO TO 7 

2 CALL RKSTEP(YN,YDN,TN,Y3H,YD3H,H*1,LEG) 

CALL RKSTEPf YN * YDN » TN , YM » YDM » H/3 • ,0,LEG) 

CALL RKSTEPI YM,YDM, TN+ ( H/3 . ) , YM , YDM , H/3 * , 0 , LEG ) 

CALL RKSTEP ( YM , YDM , TN+ ( H/3 . ) *2. , YM , YDM ,H/3. ,0,LEG ) 

DO 3 1=1,6 

EY( I ) = • 1 2 5 E-l * t YM ( 1,1) — Y3H( 1,1) ) 

3 EYD( I ) = • 12 5E- 1* ( YDM ( 1,1) -YD3H ( 1,1)) 

EV2MIN=1.E-13*(YDM( 1 *1 )**2+YDM (2,1) **2+YDM (3*1) ** 2 ) 

FVL2=AMAX1 ( FVT ,FV2MIN ) 

R = ( EYD ( 1 ) **2 + EYD ( 2 ) **2 + EYD ( 3 ) ** 2 ) /FVL2 

4 DO 5 1=1,6 

YN( I *1 )=YM( I , 1 ) +E Y ( I ) 

5 YDN ( I , 1 )=YDM( I ,1 )+EYD( I ) 

DO 6 1=1,6 

DO 6 J = 2 » LEG6 
YN ( I , J ) =Y3H( I ,J ) 

6 YDN ( I * J ) =YD3H ( I ♦ J ) 

I F ( ABS(H) .GE. ABS(TF-TN) ) GO TO 8 
TN = TN + H 

I F ( R • L T « 0,04) R=,04 
H = H /P** ,12 5 

7 IF ( ARS(H) . GT , ABS(TF-TN) ) H = T r -TN 
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I F ( A8S ( H ) .GT. HMAX) H=S IGN ( HM AX » H ) 

GO TO 2 

IF(LEG.LT.LEGMAX) GO TO 10 
CALL YDDRHS(YN,YD3H,1.»ATP(LEG+1»1) ,0) 
DO 11 1 = 1,3 

Z( I»6+LEGMAX)*YDN(I»1) 

Z< I +3 ♦ 6+LEGM AX ) =YD3H( 1,1 > 

Z< I +6 , 6+LEGM AX ) = YDN ( 1+3,1) 

11 Z( I+9,6+LEGMAX) =YD3H(I+3»1 ) 

10 DO 12 1*1,3 

XF ( I ) =YN (1,1) 

XF ( 1 + 3 ) = YDN ( I »*1 ) 

OF ( I ) = YN { 1+3,1) 

OF ( I +3 ) =YDN ( 1+3,1) 

DO 12 J=2,LEG6 
Z ( I , J — 1 ) = YN ( I ,J) 

Z( I +3 , J-l ) =YDN ( I , J) 

Z I I +6 * J— 1 )=YN{ I +3 » J~) 

12 Z( 1+9, J-l )=YDN( 1+3, J) 

AM* AM— BM AX* ( ATP ( LEG+1 ,1 )-ATP(LEG,l ) ) 
RETURN 
END 


SUBROUTINE RKSTEP ( YN , YDN , TN , YN1 , YDN 1 »H »N , L EG ) 

C THIS PROGRAM ADVANCES YN AND YDN BY A STEP OF SIZE H' TO YN1 AND YDN1 
C USING A FOURTH-ORDER RUNGE-KUTTA NUMERICAL INTEGRATION SCHEME. IF N 
C IS POSITIVE, ALL ELEMENTS OF THE MATRICES YN AND YDN ARE ADVANCED. 

C OTHERWISE ONLY THE FIRST COLUMN OF EACH MATRIX IS UPDATED. 

DIMENSION YN( 6,13) , YDN ( 6,13 ) ,D1(6,13),D2(6,13),D3(6,13)»Y(6,13) 

1 ,YN1 ( 6,13 ) , YDN 1 (6,13) 

JMAX=1 

IF(N.GT.O) JMAX=6+L£G 
H2=. 5*H 

CALL YDDRHS(YN»D1 ,H,TN,N ) 

DO 1 J= 1 , JMAX 
DO 1 1*1,6 

1 Y ( I , J ) = YN ( I ♦ J )+H2* ( YDN ( I ,J) +.25*01 ( T,J) ) 

CALL YDDRH5 ( Y »D2 »H,TN+H2 ,N) 

DO 2 J=1 , JMAX 
DO 2 1=1,6 

2 Y ( I , J ) =YN ( I , J )+H* ( YDN ( I ♦ J) + .5*D2 ( I , J ) ) 

CALL YDDRHS(Y,D3,H,TN+H,N) 

DO 3 J= 1 , JMAX 
DO 3 1*1,6 

YN1 ( I , J) =YN( I , J ) +H* ( YDN ( I , J ) + . 1 666666 7* ( D1 ( I » J ) +? • *D2 ( I , J ) ) ) 

3 YDN 1 ( I , J) =YDN ( I , J )+. 16666667* ( 01 ( I , J)+4.*D2 ( I ♦ J ) +D3 ( I , J) ) 

RETURN 

END 


SUBROUTINE YDDRHS ( Y , YDD, H , T ,N ) 

COMMON /C CP TNJ/UK »LEG»ATP( 7,2 ) , AM , 3MAX , CEX V 
DIMENSION Y(6 ,13) ,YDD(6,13) »R(3),U(3),B(6,6) 

C COMPUTE BASIC QUANTITIES COMMON TO MANY COMPONENTS OF YDD. 
LFG6=6+LFG 
DO 1 1*1,3 
R( T ) =Y( I ,1 ) 

1 U ( I ) =Y ( 1+3,1) 

AMI = AM- ( T-ATP ( LEG » 1 ) )*BMAX 
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R2 = 1 • / ( R ( 1 ) *R ( 1 ) +R ( 2 ) *R ( 2 ) +R ( 3 ) *R ( 3 ) ) 

U2 = 1./(UU)**2+U(2 )**2+U<3)**2 > 

RU=R ( 1 )*U(1)+R(2) *U ( 2 ) +R ( 3 ) *U I 3 ) 

ALPHA=-H*UK*R2*SQRT(R2) 

BETA=H*SQRT ( U2 ) *CEXV*BMAX/AM1 
GAMMA=-3.*ALPHA*R2*RU 
C COMPUTE RDD AND UDD. 

DO 2 1=1*? 

YDD ( I * 1 ) *R ( I ) *ALPHA+U( I ) *BETA 

2 YDD ( 1+3,1 )=R( I )*GAMMA+U< I ) *ALPHA 

C DEC I DF WHETHER WDD TS REQUIRFD AT THIS TIME. 

21 IF(N.LF.O) RETURN 

C COMPUTE ADDITIONAL QUANTITIES COMMON TO MANY COMPONENTS OF WDD. 
DELTA=-3.*ALPHA*R2 
EPSIL=— BETA*U2 
ZETA=-5.*GAMMA*R2 

C COMPUTE THE MATRIX B NEEDED IN THE MATRIX EQUATION WDD=B*W. 

DO 3 J= 1 * 3 
RRJ=DELT A*R ( J ) 

RUJ=EPSIL*U( J ) 

UR J = ?ET A*R ( J ) +DELTA*U( J ) 

UUJ=RRJ 
DO 3 1=1,3 
0=0 » 

IF { I .EQ. J) Q=1 . 

B( I, J)=R(I ) *RRJ+Q*ALPHA 

B( I ,J+3)=U< I )*RUJ+Q*BETA 

B( I +3 » J ) *R ( I )*URJ + U( I )*UUJ+Q*GAMMA 

3 B ( I +3 , J+3 ) =B ( I , J ! 

C PERFORM THE MATRIX MULTIPLICATION B*W TO GET WDD. 

DO 5 1=1,6 
DO 5 J=2,LEG6 
SUM= 0 , 

DO A K = 1 , 6 

4 SUM=SUM+R (I «K)*Y(K»J) 

5 YDD ( I , J5=SUM 

I F ( LEG « LE « 1 ) RETURN 
RDDMB=-3ETA*BMAX/AM1 
DO 6 J= 8 , LEG6 
DO 6 1=1,3 
Jl= J-7 

6 YDD( I »J)=YDD( I ,J) — S IGN( RDDMB.ATP ( J1 ,2 ) )*U(1) 

RETURN 

END 


SUBROUTINE BVEVAL ( XF,QF ,Z ,C , E ,DC ) 

C 

C THIS IS VERSION A OF SUBPROGRAM BVEVAL, DECK NAME ' 0PDK3 A 1 * 

C THIS VERSION DEALS WITH A FIVE-CONSTRAINT RIGHT-END BOUNDARY-VALUE 
C PROBLEM WHERE THE FIVE CONSTRAINED FUNCTIONS ARE THE THREE COMPON- 
C ENTS OF THE ORBITAL ANGULAR VELOCITY VECTOR IR CROSS V) AND THE FIRST 
C TWO COMPONENTS OF THE VECTOR WHOSE DIRECTION IS THE DIRECTION OF 
C PERICFNTER AND WHOSE MAGNITUDE IS THE ORBITAL ECCENTRICITY. THUS, IN 
C EFFECT, ALL OF THE SIX CLASSICAL ORBITAL ELEMENTS ARE CONSTRAINED 
C EXCEPT THP MEAN ANOMALY, WHICH IS FREF. 

C 

COMMON /CCPINJ/UK»LEG»ATP(7»2 ) ♦ AM , 5MAX , CEXV 
DIMENSION XF (6) , C < 1 2 ) »DC ( 1 2 ) ♦ G ( 7 ,6 ) ,Z(12,12) »E(12,13) ,R(3) *V(3) 

1 ,0F{ 6) , DUMMY! 12) 
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LEG6=6+LEG 
DO 1 1=1*3 
R ( I ) =XF ( I ) 

1 V ( I ) =XF ( I +3 ) 

R2=R (1 )**2+R(2)**2+R (3)** 2 
G( 1 » 2 ) =V( 3 ) 

G ( 1 » 3 ) =— V ( 2 ) 

G ( 2 ♦ 3 ) =V ( 1 ) 

G ( 1 * 5 ) =-R ( 3 ) 

G ( 1 ♦ 6 > =R < 2 ) 

G ( 2 * 6 ) =— R ( 1 ) 

RM=SQRT(R2 ) 

R3=RM*R2 

C1 = -1*0/RM+ CV( 1 )**2+V< 2) **2 + V< 3 )**2 ) /UK 
C2=-(R(1)*V(1)+R(2)*V(2)+R(3)*V(3) ) /UK 
DO 4 1=1,3 
DO 3 J = 1 » 3 
I F ( I .LE.J) GO TO 2 
G ( T « J ) =-G ( J» I ) 

Gt I * J+3 ) =-G ( J » I +3 ) 

2 G(I+3,J)=R(I)*R(J) /R 3-V ( I )*V( J) /UK 

3 G ( 1+3 , J+3) =(R ( I )*V( J )*2.-V<I)*R ( J) ) /UK 
G( I * I )=0. 

g ( r, r +3 ) =o * 

G( 1+3 , T )=G( 1+3 , I ) +C1 

4 G( 1+3 , 1+3 ) =G( 1+3 , 1+3 )+C2 
DO 5 1=1,3 

DC ( I )=0. 

DO 5 J = 1 , 3 

5 DC( I) =DC( I }+G ( I , J) *R ( J) 

DO 8 1=1*2 
SUM=0. 

DO 7 J=1 , 3 

7 SUM=SUM+G(I,J)*DC( J) 

8 DC( 1+3 )=C( 1+3 )+R( I ) /RM+SUM/UK 
DO 9 T = 1 ,3 

9 DC( I ) =C( I )-OC( I) 

DO 10 1=1,5 
DO 10 J=1,LEG6 
E ( I , J)=0, 

DO 10 K = 1 » 6 

10 E(T,J)=E(I,J)+G(I » K ) *Z ( K » J ) 

R 5 = X F ( 1 ) *QF( 1 )+XF ( 2 ) *QF ( 2 ) +XF ( 3 ) *QF ( 3 ) 

C3=-UK/R3 
C4=-3. 0*C3*RS/R2 

DC ( 6 ) = XF(4)*QF(4)+XF(5)*QF(5)+XF(6) *OF ( 6 ) -C3+RS 

DO 16 1=1,3 

DUMMY ( I ) =— C 3 *QF ( I )-C4*XF( I ) 

DUMMY ( 1+3 )=QF( 1+3 ) 

DUMMY ( 1 + 6 )=-C3*XF( r ) 

16 DUMMY ( 1+9 ) =XF ( 1+3) 

DO 17 J=1,LEG6 

DO 17 K = 1 ,12 

17 E( 6 ,J) =-DUMMY ( K) *Z(K*J) + E(6,J) 

RETURN 

END 


SUBROUTINE AD JUST ( 00 , E ♦ Z , DRF , DVF , CK , LEGMAX , ATP ) 
COMMON/ WADJ/A( 8) ,DXF(6) 
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DIMENSION Q0(6) » E ( 1 2 » 1 3 ) » Z { 1 2 * 1 2 ) , ATP (7, 2) 

CALL SOLVE(E*LEGMAX) 

L6=6+LEGMAX 
L7=LEGMAX+7 
DO 5 1=1*6 
DXF ( I > =0 . 

DO 5 K=1»L6 

5 DXF { I )=DXF( I )+Z ( I »K)*E(K»7+LEGMAX) 

DRF=SORT ( DXF ( 1 ) **2+DXF (2 ) **2+DXF(3 ) **2 ) 
DVF=SQRT(DXF(4)**2+DXF(5 ) **2+DXF ( 6 ) **2 ) 

DU2 = S0RT(E( 1 , 7+LEGMAX )**2+E< 2 » 7+LEGMAX ) **2+E < 3 .7+LEGMAX )**2) 
DUD2=SQRT(E(4»7+LEGMAX)**2+E< 5 , 7+LEGMAX ) **2 + E ( 6 » 7+LEGMAX ) **2 ) 
A ( 1 ) = • 2/DU2 
A(2)=.0003/DUD2 
A ( 3 ) = 1 • 0 

IF(ATP(l,2).LT.O) A ( 3 ) = ( . 5* ( ATP ( 2 * 1 ) -ATP ( 1 , 1 ) ) /ABS ( E ( 7 * L7 ) ) ) 
CK = AM I N 1 ( 1 • 0 » A ( 1 ) * A ( 2 ) * A ( 3 ) ) 

DO ,8 1=2 » LEGMAX 
12=1+2 

A( I2) = ( . 5* ( A TP ( 1 + 1*1) — A T P ( 1*1) ) / A RS ( E ( I+6»L7)-E(I+5*L7) ) ) 

8 CK=AMIN1(CK*A( 12) ) 

DO 6 1=1,6 

ATP ( 1+1,1 )=A TP ( I+1,1)+CK*E (I+6»7+LEGMAX) 

6 Q0( n=QO(T)+CK*E(I, 7+LEGMAX) 

UM=SQRT (00(1)** 2+00(2) **2+Q0( 3 ) **2 ) 

DO 7 1=1*6 

7 Q0( I )=00( I)/UM 

RETURN 
END 


SUBROUTINE SOLVE ( A , LEGMA X ) 

DIMENSION A ( 12*13) 

L6=6+LEGMAX 
L 7 = 7+L EGMAX 
DO 6 N = 1 * L6 
I B I G=N 
DO 1 I =N * L6 

I F ( AB S ( A ( I , N ) >.GT. AB S ( A ( I B I G , N ) ) ) IBIG=I 

1 CONTINUE 

IF ( IBIG.EQ.N ) GO TO 3 
DO 2 J = N * L7 
0=A ( N, J) 

A ( N * J ) =A ( IBIG.J) 

2 A ( I B I G , J ) =0 
3 DO 5 1 = 1, L6 

I F ( I.EO.N) GO TO 5 
0= A ( I ,N)/A(N,N) 

M = N+ 1 

DO 4 K = M ♦ L7 

4 A ( I , K ) =A ( I * K ) -Q*A ( N ♦ K ) 

5 CONTINUE 

6 CONTINUE 

DO 7 1=1, L6 

7 A( I , L7 ) = A ( I » L7 ) / A ( I ,1 ) 

RETURN 

END 


DATA FOR CASE 1 
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1 

398601.56 12644651. 22384.40794.1540729 

5087.493 4132.6342 114.70552 4.900796 

.4668414 -.7801086 -.1528746 -.5666 E- 

1674.843496091.504 -87016.17 0.0 
450 

1034. 1.0 

1434. -1.0 

1689.8176 1.0 

20417.2732-1.0 
20546.38321.0 


100 . 

-6.0577962.0669144 
-.6626 E-3-,1548 E-3 
0.0 
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CASE NUMBER 1 


UKa 390601.56 INITIAL MASSa T . 264465 1 OE+07 MASS RATE* 2.23844079E*04 EXHAUST VELOCITYa 4. 15407290E*00 HMAXa 100.000 


IMAXx 50 LEGMAX- 4 



TIME 

TYPE OF ARC ( * COAST.- BURN) 




1 . 03400000E*03 

1.00000000E*00 




1 .434000002*03 

. . -1*000000002*00 - — 




1 .68981 76OE*03 

1.00000000F*00 




2.0417273?E*04 

-1.00000000E*00 




2.0546383?E*04 

1.00000000E*00 



INITIAL STATE XO 

5087.4930fi00 

4132.6342000 114.7055200 4.9007960 

-6.0S77962 

.0669144 

ESTIMATED QO 

.4668414 

... -.7801086- - - -*152GZ46 — **0005666 

..*0006626 

-*0001548 

DESIRED FINAL C -1 

,674A4340E*03 9. 

60915040E*04 -8.70161700E*04 0. 0. 

-0. 

-0. 

ITERATION NUMBER 

1 





COAST ABC 

LEG- 1 STATE AT END 6.408976E*03 1.340597E+03 1.277190E.02 1.5796512*00 -7.638430E*00 -3.liS99SE.03 

COSTATE AT END 1.678606E-01; - -9.661 125E-6* -1*961 1 lSE-01 -8.810358E-O4 -2.277O36€-04~ -5»677Gi4fr-05 

PSYa 6.104706E-02 ALPHA- -6.08B936E*01 CALCULATED COAST TIME- 4.0000002*02 

SEMIMAJOR AXISa 6.546326E*03 RMINa 6.53ll?0E*03 RMAX- 6.56l4B2E*03 ENERGY- -3.0444682*01 
PEBIOO. 5.271167E*03-HMAOa — 5*10019 1 F*04- -H -VECTOR 9*713954E*0? 2.2172162*02 -5.1072192*04- 

EMASa 2.315202E-03 E VECTOR 7.232001E-05 -2.314056E-03 -8.670576E-06 RNAG- 6.548930E+03 
BURN ARC 

LEG- 2 STATE AT END 6.535795E+02 -9.242956E*02 4*411060E*01 -6.737298E-01 -1.018129E*01 -$*422B21C-«1 

COSTATE AT END -6.29U97E-02 -9.793821E-01 -2.014162E-01 -B.896451E-04 1.235212E-04 1.567432E-0S 

MASS AT END OF LEG- 6,9183252*06 

COAST AOC . ■ - 

LEG- 3 STATE AT END -4.213018E*04 -1 .570124E*03 -8.0565l7E*02 -7.446084E-02 1.591464E*00 8.16T224E-0* 

COSTATE AT END -2.683806E-02 2.094087E-01 9.794587E-01 2.54i943E-05 1.117771E-06 T.909943E-0T 

PSY- 7.576405E-01 ALPHA- -1.636083E*0l CALCULATED COAST TIME- 1.872746E*04 

SEMIMAJOR AXIS- 2.436317E«04 RM1N- 6.558651E*03 RMAX- 4.216769E«04 ENERGY- -8. 18Q413E*00 
PERIOD- 3.784522E+04 HMAG- 6.726664E*04 H VECTOR 1.153930E*03 3.500656E*03 -6.716557E*04 

EMAG- 7.3O7965E-01 E VECTOR 7.302401E-01- 2.492529E-02 1.3B4498E-02 RMAO- 4.216T13E*04 

BURN ARC 

LEG- 4 STATE AT END -4.214131E*04 - 1 . 336984E*03 -6.660099E*02 -1 . 017078E-01 2.062185E*00 2*2T6494E*00 

COSTATE AT END -2.35S503E-02 2.095438E-01 9.795174E-01 2.543611E-05 9.739700E-07 1.166818E-0T 

MASS AT END OF LEG- 4,028275E*06 

SEMIMAJOR AXIS- 4.217548E+04 RMIN- 4.216691E*04 RMAX- 4.218406E.04 ENERGY- -4.725513C*00 
PERIOD- 8.619865E*04 HMAG- 1.296S81E*05 H VECTOR -1.672874E*03 9.608644E*04 -8.7039132*04 

EMAG- 2.032591E-04 E VECTOR -1 .798668E-04 -6.525732E-05 -6.858351E-05 RMAGa 4.216777E*04 
TOTAL BURN TIME- 3.849276E+02 ARC TIMES 4.000000E*02 2.558176E*0? 1.872746E*04 1.291100E*02 

DC -1 ,96891 1E+00 5.063665E*00 2.296253E*01 1.798668E-04 6.525732E-05 1.145031E-08 

9.004042E-07 1.874982E-08 -9. 780427E-05 0. 0. 0. 

DETERMINANT OF E- -4.953781E+06 OlAGONAL OF E 2.363277E*06 7.716025E*04 2.779820E*04 5.059089E*03 

-1.2Bp?83E*03 1 . 344838E*0 1 1.453115E-02 -9.844767E-06 5.183974E-03 -1 .510309E-08 

OXF -1.395360E-01 -2.832318E-01 3.898002F-01 2.955619E-04 -5.430284E-04 1.182267E-04 
DRF- 5.016321E-01 DVF- 6.294555E-04 CK- 1.000000E*00 EVT- 1.000000E-14 
CK-MIN OF 9.25126E+02 1.46230E*02 1.00000E*00 7.52355e*05 4.49457E*03 2.23327E*05 

CHANGE REQUISTED IN INITIAL COSTATE. SWITCHING TIMES AND FINAL TIME 


A21 


NEW 00 5.0M979E-01 -8.463290E-01 -1 .650O37E-O1 -6. 1 44147E-04 -7. 187790F-04 -1 .680213E-04 

NEW SWITCH TIMES 1.433S32E+03 1.689650E+03 2.041919E*04 2.0548305*04 

ITERATION NUMBER 2 

COAST ABC 

LEG* 1 STATE AT END 6.408711E+03 1.34188lE*03 1.277195E+02 1.581180E*00 -7.6381 i 0E*00 -3.085522E-03 

COSTATE AT ENO 1.820784E-01 -1 . 048089E*00 -2. 127385E-01 -9.556442E-04 -2.473559E-04 -6. 175324E-05 

PSY* 6.10213BE-02 ALPHA* -6.088936E*01 CALCULATED COAST TIME* 3.998319E*02 

SEMIMAJOR AXIS* 6.546326E*03 RMIN* 6.531170E*03 RMAX* 6.561482E*03 ENERGY* -3.044468E*01 
PERIOD* 5.271 167E+03 HMAG* 5. 108191E+04 H VECTOR 9.713954E*02 2.217218E+02 -5. 107219E*04 
EMAG* 2.315202E-03 E VECTOR 7.232001E-05 -2.314056E-03 -8.670576E-06 RMAG* 6.548933E*03 
BURN ARC 

LEG* 2 STATE AT END 6.535921E*03 -9,229923E*02 6.411993E*01 -6.722314E-01 -1.01814SE*01 -5.422219E-01 

COSTATE AT ENO -6.826R09E-02 .062571E*00 -2. 185335E-61 -9.652236E-04 1.336812E-04 1.684696E-05 

MASS at END of leg* 6.918322E*06 

COAST ARC 

LEG* 3 STATE AT END -4.213o35E*04 -1 .570370E*03 -8.0527l2E*02 -7.415446E-02 1.5914491*00 8.i47749E-02 

COSTATE AT ENO -2.91 1961E-02 2.269079E-01 1.062597E*00 2.757984E-0S 1.191153E-06 8.667272E-07 

PSY* 7.576990E-Q1 ALPHA* -1 ,6360T8E*01 CALCULATED COAST TIME* 1.8729S4E404 

SEMIMAJOR AXIS*. 2.436324E*04 RMIN*. -4.5S86S1C*03 RMAX* 4.2147«3E*04 ENERGY* -8.1803881*00 
PERIOD* 3.784539E*04 HMAG* 6.726666E*04 H VECTOR 1.153300E*03 3.500816E*03 -6.716560E*0i 
EMAG* 7.307973E-01 E VECTOR 7.302390E-01 2.49B246E-02 1.384107E-02 RMAG* 4.216730E*04 

BURN ARC 

LEG* 4 STATE AT END -4.214144E*04 -1 .33726lE*03 -6.65621 1E*02 -1.014U7E-01 2.061643E*00 2.278624E*00 

COSTATE AT END -2.555755E-02 2.270517E-01 1.062662E*00 2.759790E-05 1.035336E-06 1.373526E-07 

MASS AT ENO- OF LEG - 4 t 028264P* 0 6 

SEMIMAJOR AXIS* 4.216B17E«04 RMIN* 4.216790E*04 RMAX* 4.216844E«04 ENERGY* -4.726332E*00 
PERIOD* 8.617623E*04 HMAG* 1.296468E*05 H VECTOR -1 .6748425*03 9.609199E*04 -8.70i623E*O4 

EMAG* 6.319555E-06 E VECTOR -6.31S798E-06 -9.661769E-0B 1.492601E-08 RMAG* 4*214790E*04 

TOTAL BURN TIME* 3.849281E*02 ARC TIMES 3.998319E*02 2.558l78E*02 1.872954E*04 1.291103E*02 

DC -l.123247E-03._-4.62563T.E-Ql . 5. 628240 E -02 6J18748E-06 9.661769E-08 -2.406201E-11 

-6.245005E-17 4.207506E-11 1.198258E-05 0. 0. 0. 

DETERMINANT OF E* -4 ,952307E*06 DIAGONAL OF E 2.178922E*06 7.jl909lE*04 2.557792E*04 4.666866E*03 

-1.39o846E*03 1.243308E*01 1.4541S9E-02 -1.16738SE-0S 5.562259E-03 -1 .438006E-08 

DXF -8.877080E-03 -8.556068E-03 -1 .748166E-03 -7.214129E-07 -1 .813445E-06 -1.194664E-05 
ORF* 1 .245251E-02 DVF* J.210501E-05 CK* 1.000000E*00 EVT* l.OOOOOOE-14 
CK*MIN OF 1.06&15E+05 1.604975*04 1.00000C*00 l.S0719E*04 2*52071E*07 1.10879E*05 

change REQUISTED IN INITIAL COSTATE. SWITCHING TIMES AND FINAL TIME 

1.556672E-07 4.476341E-07 -1 .81 1483E-06-4.608971E-10 1.774056E-10 -l .587006E-09 

-2.636723E-05 4.441062E-05 4.159234E-04 -1 .662874E-04 

END OF ITERATION NUMBER 2 

NEW QO 5.061981E-01 -8.463286E-01 -1 .65805SE-01 -6.144152E-04 -7. 187788E-04 -1 .680229E-04 
NEW SWITCH TIMES 1.433832E+03 1.689650E*03 2.041919E*04 2.054830E*04 



SUMMARY tables 

ITERATION total burn time length OF burn and coast arcs 

number 

1 3.8492760E*02 4.0000000E+02 2.5581T60E*02 1 ,8727456E*04 1 .2911000E.02 

2 3.8492806E*02 3.9983188E*02 2.55817776*02 1.87295396*04 1 *291 10296*02 

iteration ERROR IN BOUNDARY C0N0ITI0NS-DCU-8) 

number 

1 -1.968911E*00 5.063665E*00 2.296253E*0l 1.798668E-04 6.525732E-05 1.1456316-08 9.0040426-07 1.8749826-08 

2 -1.123247E-03 -4.829637E-01 5.628740E-02 6.3187986*06 9.6617696-08 -2.406201E-1 1 -6.245005E-17 4.2075046-11 

ITERATION " DC (9-12) DRF 0V6 CK 

NUMBER 

1 •9.78Q4266E-05 0. 0. 0. 5.0163209E-01 6.29455S2E-04 1.00000006*00 

2 1.19825836-05 tU. 0. 0.- -U24S2f0TE-O2 W21050116-05 1.00000006*00 



